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,—1 ■ ABSTRACT 

: 

• Context. In the last couple of years a population of very massive (M* > 10^^ Mq), high-redshift (z > 2) galaxies has been 
^ ' identified, but its role in galaxy evolution has not yet been fully understood. 
• i-H , Aims. It is necessary to perform a systematic study of high-redshift massive galaxies, in order to determine the shape of the 
' very massive tail of the stellar mass function and determine the epoch of their assembly. 
%-i _ Methods. We selected high-z massive galaxies at 5.8pim, in the SWIRE ELAIS-Sl field (1 deg'^). Galaxies with the 1.6pim stellar 
. . . . peak redshifted into the IRAC bands ~ 1 — 3, called "IR-peakers" ) were identified. Stellar masses were derived by means of 
spectro-photometric fitting and used to compute the stellar mass function (MF) at z = 1 — 2 and 2 — 3. A parametric fit to the 
MF was performed, based on a Bayesian formalism, and the stellar mass density of massive galaxies above z = 2 determined. 
Results. We present the first systematic study of the very-massive tail of the galaxy stellar mass function at high redshift. A 
total of 326 sources were selected. The majority of these galaxies have stellar masses in excess of lO'^^ Mq and lie at z > 1.5. 
The availability of mid-IR data turned out to be a valuable tool to constrain the contribution of young stars to galaxy 
SEDs, and thus their M^/L ratio. The influence of near-IR data and of the chosen stellar library on the SED fltting are also 
discussed. The z = 2 — 3 stellar mass function between 10^^ and ~ 10^^ Mq is probed with unprecedented detail. A signiflcant 
evolution is found not only for galaxies with M ^ 10^^ M©, but also in the highest mass bins considered. The comoving 
number density of these galaxies was lower by more than a factor of 10 at z = 2 — 3, with respect to the local estimate. 
SWIRE 5.8/im peakers more massive than 1.6 x 10^^ Mq provide 30—50% of the total stellar mass density in galaxies at z = 2 — 3. 



Key words. Galaxies: evolution - Galaxies: mass function - Galaxies high-redshift - Galaxies: fundamental parameters - Galaxies: 
statistics - Infrared: galaxies 



1. Introduction 
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ogy, particularly controversial when dealing with massive 
(^stars > 10" Mo) objects. 

The assembly of massive galaxies is one of the crit- 
ical questions in the cosmic evolutionary scenario. The 
uniform properties of local early-type galaxies and of the 
fundamental plane have inspired the so called "monolithic 
collap se" scenario ((Eggen et al.l Il962 ; Chiosi &: Carrarol 
20021 ) , in which galaxies formed in the remote past through 
huge events of star formation and subsequently evolved 
passively across cosmic time. On the o ther hand, in 
the more recent "hierarchical" scenario ( White fc Rees 
1978": 'Kauffmann et al."l996 HKauffmann fc Charlotlll998t 
Somcrvillc & Primack 19991), massive galaxies assemble by 
mergers of lower-mass units, with the most massive ob- 
jects being born in the latest stages of evolution, at z < 1. 

The availability of several powerful tracers of star for- 
mation (e.g. UV continuum, optical recombination lines, 
far-IR emission, sub-mm light) has favored the popu- 
larity of studies of the comoving st ar formation density 
( Madau et al.lll996l : Lillv et aTlll996[ ) in galaxies at various 
redshifts. It is now well determined that the Universe ex- 
perienced an epoch of enhanced star formation in the past, 
peaking at z ~ 1 — 2, with a subsequent decline o f at least 



one ord er of magnitude to the present time fe.g. iHopkins 



2004; Rudnick et aIll2003HFiores et al.|[l999HMadau etal 
19961 among others) 



An alternative approach consists of studying the mass 
already assembled in galaxies, instead of the amount of 
stars being formed. The integral of the past star formation 
density provides the stellar mass density at a given epoch: 
a complementary constraint on cosmic galaxy evolution. 

Thus, the build up of the stellar mass across cosmic 
time has become one of the major topics in observa- 
tional cosmology, and has overtaken the classic Madau- 
Lilly diagram as the central tool for studying galaxy evo- 
lution. The large observational effort dedicated to this 
subject has shown that the global stellar mass density 
incre ases from early epochs to t he low-redshift Universe 
e.g. iBrinchmann fc E lii3 I2 OO0I: lOickinson et al.l l2003a : 
Fontana et a .120031 12004. 2006[ lRudnick et al.ll2003l . l2006 ; 



Drorv et all |2005| ). Very deep surveys have been ex- 
ploited to describ e the shape of the stellar mass function 
at high redshift (iFontana et all [200I iDrorv et all 120051: 
Gwvn fc Hartwickll2005( r" but a clear picture on the role 
of very massive galaxies has not yet emerged. 

Several pieces of evidence exist that fully formed mas- 
sive galaxies were already in place at redshift z ~ 2 — 3. 

A substantial population of luminous red galax- 
ies at redshifts z > 2 (known as "distant red 
galaxies", DRGs) was found in the Faint InfraRed 



Extragalactic Survey CFIR ES. iFranx et all l2003h. Based 



on near-IR spectroscopy (|Forster Schreiber et al. 2004 : 



van Dokkum et al. I2OO 



these galaxies turned out to 
be massive (M = 1 — 5 x 10-^^ M©), evolved (ages of 
1 — 2.5 Gyr) systems, probably descendants of galaxies 
whic h started forming at r edshift z > 4. Based on FIRES 
data, Rudnick et al. ( 20031 ) inferred that DRGs contribute 
^ 50% of the global stellar mass density at z = 2 — 3. 



Dickinson et al.l (j2003a|) exploited Hubble Deep Field 
North NICMOS data to derive the stellar masses of galax- 
ies up to z = 3. The study of the global stellar mass den- 
sity highlighted that 50 — 75% of the present-day stellar 
mass was already in place at z 1, while only 3 — 14% 
had been already assembled at z ~ 2.7. 

Direct determinations of the galaxy mass func- 
tion based on near-IR deep imaging by the Spitzer 
Space Telescope indicate that the number of massive 
galaxies does not significantly e v olve up to at least 



1 dFranceschini et al. 20061 : Fontana et al. 2006t 



Bundv et alJl2005^ . 



Usin g data from 



(GDDS, [Abraham et al 



the Ger nini Deep Deep Surve; 
20041) . iGlazebrook et all ()200 



identified a population of red {[I — K] > 4 Vega mag.) 
galaxies at z ~ 2 with stellar masses in excess of lO""^""^ 
Mq. These objects contribute roughly 30% of the to- 
tal stellar mass density of the Universe at that epoch. 
McCarthv et alJ ( 20041 ) estimated the age of red galaxies 



at z = 1.3 — 2.2 in the GDDS, deriving a median age 
of 1—3 Gyr, and a star formation history dominated by 
very powerful bursts (300-500 Mq yr'^). These massive 
galaxies must have undergone a rapid formation process 
at z > 1. 



Exploit ing data from the K20 survey (jCimatti et al 



2002al lblO). rDaddi et all (|2004l ) identified few luminous K 



band selected galaxies atl.7<z<2.3 with stellar masses 
M ~ 10" - 5 X 10" [Mq]. C ombining deep K20 s pec- 
troscopy and HST-ACS imaging. ICimatti et all (j2004l ) dis- 
covered four old, fully assembled spheroidal galaxies at 
1.6 < z < 1.9: the most distant such objects currently 
known. The stellar mass of t hese galaxies turned out to 
be in the range 1-3 x 10" Mg. lFontana et alJ (|2004l ) stud- 
ied K20 galaxies as well, showing that massive (M > 10"'^^ 
M0) galaxies are easily found up to z ~ 2. These authors 
also report on the stellar mass function: only mild evolu- 
tion (~2— 30%) is detected to z = 1, but only ^ 35% of 
the z = stellar mass locked up in massive objects was 
assembled by z = 2. 

At even higher redshifts, Rigopoulou et al. ( 20061 ) have 
recently studied a population of z ~ 3 Lyman-break galax- 
ies (LBGs) with st ellar masses in excess of 10^^ M©. 
McLure et al. I (l2006^ identified nine LBGs at z > 5 in the 



UKIDSS ( Lawrence et al.l[2006l ) survey, over an area of 0.6 
deg^ . A stacking analysis suggests that t he typical stellar 



Mobasher et al 



mass of these sources is > 5 x 10^° Mq 

( 20051 ) analyzed the proper ties of J-dropouts in th e Hubble 
Ultra Deep Field (HUDF, iBeckwith et"alll2006f) . exploit- 
ing Spitzer photometry. They identified a z ~ 6.5 candi- 
date that was interpreted as a post-starburst galaxy with 
a surprisingly high ste llar in ass of 5.7 x 10^^ Mq (but see, 
for example, Yan et al . l2004l for a different interpretation). 



Drorv et all ( 20051 ) studied the stellar mass function of 



galaxies in the FORS Deep Field (FDF iHeidt et alJ l2003h 
and GOODS/CDFS (|Giavahsco et al.ll2004l ) fiekTover a 
total area of 90 arcmin^ and found that the total stellar 
mass density at z = 1 is 50% of the local value. At z = 2, 
25% of the local mass density was already assembled, and 
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at z — 3 and z = 5, at least 15% and 5% of stellar mass, 
respectively, was already in place. Massive (Af > 10^"'^ 
Mq) galaxies existed over the whole redshift range probed, 
up to z = 5. The number density of these massive galaxies 
evolves very similarly to galaxies with Af > 10^" M©, 
decreasing by 0.4 dex to z = 1, 0.6 dex to z = 2, and 1 
dex to z = 4. 

By analyzing the properties of if -selected galaxies 
in the ^ 131 arc min^ of the GOODS-CDFS survey, 
Caputi et al. I (l2006h found that the vast majority (85 — 
90%) of local M > 2.5 x 10" Mq galaxies appears to be 
already in place at z ~ 1. These authors also infer that 
roughly 65 — 70% of these galaxies assembled at z = 1 — 3 
by means of obscured, intense bursts of star formation, 
while the remaining could be in place at even higher red- 
shifts (z = 3-4). 

The observational challenge is that large volumes are 
needed to find representative samples of such rare very 
massive galaxies at high redshift. The Spi tzer Wide-area 



InfraRed E xtragalactic survey (SWIRE, iLonsdale et al 
20031 12004 observed ~49 deg^ in the seven Spitzer chan- 
nels, and is therefore ideal to find rare objects. Its vol- 
ume is large enough to detect ~85 DM h aloes of mass 
> in the 2 < z < 3 redshift range (jjenkins et al 



2OOII: IMo fc White! [2003 ). These haloes are predicted to 
host the most luminous {Ltoi > 10^^ Lq) and most mas- 
sive {M^, > several 10^"^ Mq) galaxies ever to exist. 

The Infrared Array Camera (IRAQ . iFazio et al.ll2004h . 
onboard Spitzer ( Werner et al.l 2004), samples the rest- 
frame near-IR light of distant galaxies. A near-IR selection 
not only directly probes the low-mass stars dominating 
the baryonic mass of a galaxy, but also is minimally af- 
fected by dust extinction. Therefore Spitzer is best suited 
to the study of the stellar content of galaxies up to z = 3. 
Moreover Spitzer allows detection of galaxies that would 
be missed by rest frame UV selectio n, for example distant 
red galaxies fe.g. lDaddi et al.ll2004j ). 

We take advantage of the shape of near-IR spectral 
energy distribution (SEDs) of galaxies to identify high- 
redshift objects on the basis of IRAC colors. Our selec- 
tion is base d on the detection of the l.diim stellar pea k 
in galaxies ( Sawicki 20021 : Simpson &: Eisenhardll 1999t ). 
redshifted to the IRAC domain. It is also worth noting 
that galaxies selected in this way benefit from a negative 
k-correction in the IRAC bands, because the slope of a 
galaxy's SED is negative redward of the 1.6/Ltm peak. In 
this way we performed a systematic search for M > 10^^ 
Mq galaxies at z > 1. 

The analysis is carried out in the central square de- 
gree of the ELAIS-Sl SWIRE field, where optical, near-IR 
(J, Kg) phot ometry and optical spectroscopy are available 
(Berta et al., 20061 Dias et al., in prep.. La Franca et al., 
in prep.). This area, and the sampled volume, are big- 
ger than any other previously explored for studying very 
massive galaxies at z = 1 — 3, wh ich were limited to very 
deep, pencil-beam surveys (e.g. Dickinso n et al.l 2003at 



2006h . 



This paper is structured as follows. Section 2 presents 
the data available in the ELAIS-Sl field; in Sect. 3 we 
present our selection criterion; then Sect. 4 discusses the 
photometric estimate of redshifts. Section 5 deals with 
the estimate of the stellar mass in galaxies, and presents a 
very detailed analysis of the influence of mid-IR and near- 
IR constraints on it. Experiments with different stellar li- 
braries and IMFs are also discussed. Section 6 presents the 
stellar mass function of our galaxies, including complete- 
ness correction and a parametric fit based on a Bayesian 
formalism. Finally, Sects. 7 and 8 discuss results and draw 
our conclusions. 

Throughout this work, we adopt a standard Hq = 71 
[km s^^ Mpc^^], n„i — 0.27, JIa — 0.73 cosmology, unless 
otherwise stated. 



2. Available Data 

The ELAIS-Sl field 
-44°00'00", J2000.0) 



(RA = 00''38'"30% Dec = 
represents t he minimum of the 

1998h 



Milky Way 100 fj,m. cirrus emission (jSchlegel et al 



in the southern hemisphere, with an average emissivity of 
0.38 MJy/sterad. 



2.1. Spitzer SWIRE data 

This area was targe ted by the Spitzer Space Telescope 
(jWerner et al.ll2004h . as part of the Spi tzer Wide area 



Infra-Red Extragalactic survey (SWIRE, ILonsdale et al 



20031 l2004f). A total of ^ 7 deg^ were observed with the 



IRAC (|Fazio et al.l |2004| ) and MIPS (|Rieke et all \20Q4) 



cameras. 



Data processing is described bv ISurace et al. (2004), 
Surace et al. (in prep.), Shupe et al. (in prep.) and Afonso 
Luis et al. (in prep.). It consists of Basic Calibrated Data 
(BCD) by the SSC pipeline plus post-processing aimed 
at artifact removal, mosaicking and source extraction. 
Mosaicking was performed with the SS C routine MOPEX, 
and source extraction with SExtractor ( Bertin fc Arnouta 
Il996l ). For unresolved sources (i.e. the case examined 
here), IRAC fluxes were extracted through a 1.9" diameter 
aperture and corrected to total fluxes following SSC pre- 
scriptions; MIPS fluxes were extracted by means of PRE 
fitting (see Surace et al., and MIPS Data Handbook 2006). 

The resulting 5a depths are 4.1, 8.5, 43, 48, 400 /iJy at 
3.6, 4.5, 5.8, 8.0, and 24 ^m respectively. Observations at 
70 and 160 /um lead to 26 and 166 mJy 5cr limits (Surace 
et al., Shupe et al., Afonso Luis et al., in prep.). 

For more details on Spitzer data reduction, calibration 
and ca talog extraction see the SWIRE delivery documen- 
tation ( Surace et al. 2004 . and following release^). 



Drorv et al]|2005t ICwvn fc Hartwick,2005, : iFontana et al 



^ available at the Spitzer Science Center Legacy Program 
web page, |http://ssc.spitzer. caltech.edu/legacy/, 
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2.2. Optical ESIS imaging 

The ELAIS-Sl area is the target of extensive optical fol- 
low up carried out with ESO telescop es: the ESO-Spitze r 
Imaging extragalactic Survey (ESIS, iBerta et al. 20061 ). 
The optical ancillary data cover ~ 5 deg^ in the B, V, 
R ba nds, o bserved with the Wide Field Imager (Baade 
et al. Il999l . on the 2.2m ESO-MPI telescope in La Silla) 
and i n the /, z bands, obtained with VIMOS (Le Fevre 
et al. l2002i on the VLT). The B, V, R observations and 
re sults over the c entral 1.5 deg^ area were presented in 
3rta et al.l (|2006l) . while the VIMOS data are in an ad- 
vanced reduction stage (Berta et al., in prep.). 

The reduction of B, V, R data was performed by us- 
ing IRAF standard tasks and self-built routines. Sky flat- 
field frames were acquired during each night and applied 
to images obtained on the same date. In order to obtain 
a uniform background and photometric zeropoint across 
the field, super-sky-flat frames were built and applied for 
each observing night, taking care of masking very bright 
sources. 

Astrometric calibr ation was performed on observations 
of IStone et al. ( 1999t ) astrometric fields. In order to ac- 
curately map coordinates, a TNX algorithm was chosen, 
combining a gnomonic projection and non-linear polyno- 
mial distortions. The r.m.s coordinate difference between 
the ESIS and the GSC 2.2 catalogs turned out to be ~ 0.1 
arcsec in RA and Dec. The relative astrometric accuracy 
between different ESIS pointings is better than 0.05 arcsec 
r.m.s. 

Photometr ic calibration was obtained by observing 
Landoltl (jl992r ) standard fields, in order to calibrate pho- 
tometric zeropoints. Since observations were spread across 
several years, particular care has been taken into account- 
ing for zeropoint differences between the numerous science 
frames. Moreover, BVR color-curves were built, and mag- 
nitudes were transformed to the Johnson- Cousins stan- 
dard photometr ic system. Cat alogs were extracted by us- 
ing SExtractor (jBertin fc Arnou ts 1996). Kron total mag- 
nitudes are adopted for extended sources, while aperture 
magnitudes, corrected using the observed PSF, are used 
for unresolved objects. The final catalog reaches 95% com- 
pleteness at B V 25 and R 2 4.5. 

We defer to lBerta et al. 1 (l200d) for a further description 
of these data and their analysis. 



2.3. Near-IR J and Kg imaging 

The square degree at the center of the ELAIS-Sl area 
includes J and K,, observations carried out with the SOFI 



(jMoorwood et al.lll998f ) camera on the NTT, during four 
different periods in 2002 and 2003 (Dias et al., in prep.). 

Pre-reduction, sky-subtraction and mosaicking were 
carried out in the standard way, using the IRAF environ- 
ment. The astrometric ma pping was calibra ted by using 



aperture photometry of point-like objects. The average J 
and Ks magnitude difference for the sources in common 
with the 2MASS catalog is ^ 0.02 mag in both bands. 

Catalog extraction w as performed using SExtractor 
(jBertin fc Arnoutd Il996h . In this work we make use of 
total magnitudes. The resulting catalog is 95% complete 
to J - 19.8 and Ks = 18.73 (Vega). 

The optical, near-IR and SWIRE/Spitzer catalogs 
were matched with a simple closest-neighbo r algorithm, 
adop ting a 1.0 arcsec matching radius (see iBerta et"al 
20061) . 



the ESIS R band images (jBerta et al. I l2006h . reaching a 



0.22 and 0.16 r.m.s. uncertainty for RA and Dec, re- 
spectively. Photometric calibration was computed using 



2.4. Other imaging data 

The central area benefit s from X-ray ob servation by the 
XMM-Newton telescope (jPuccetti et al.| [2006'). The whole 
ELAIS-Sl field wa s observed at 1.4 GHz w ith the ATCA 
radio telescope bv ICruppioni et all ( 1999[ ) down to a 80 
^Jy r.m.s. and by Middelberg et al. (submitted) down to 
30 /iJy r.m.s. Finall y, the ultraviolet G alaxy Evolution 
Explorer (GALEX, iMartin et all |2005| ) Deep Imaging 
Survey (DIS) included the ELAIS-Sl field. 

2.5. Available spectroscopic data 

The ELAIS-Sl area has been spectroscopically surveyed 
by different projects, in particular focused on the central 
field 

iLa Franca et al. (2004) performed optical spec- 
troscopy of ISOCAM 15/im counterparts, with the 
2dF/AAT, ESO-Danish 1.5m, ESO 3.6m and NTT tele- 
scopes, over the spectral range between 4000-9000 A. 

During 2004 and 2005, the XMM-NIR area was the tar- 
get of 5000—9500 A low resolution spectroscopy, carried 
out with VIMOS- VLT. The primary targets of this survey 
are X-ray sources, if-selected galaxies {Kg < 18.5), and 
24 /xm SWIRE objects. Further spectroscopy of optically 
bright (i? < 21) X-ray and mid-IR sources was obtained 
at the 3.6m/ESO telescope, during FaU 2005. These ob- 
servations will be presented in future works by La Franca 
and collaborators. 

Overall, at the present time, ~1250 spectroscopic red- 
shifts are available in the ELAIS-Sl SWIRE area. 

3. Selection criteria 

In order to select high redshift galaxies dominated by stel- 
lar emission in the restframe near-IR, we have exploited 
the "IR-pe ak" technique, des cribed in Lonsdale et al. (in 
prep. ) andEertaeian (|2007l ). 

The near-IR emission of a galaxy is characterized by a 
peak centered at 1.6/xm, due to the combination of the 
Planck spectral peak of low-mass stars (dominated by 
type M), a minimum in the H" opacity in stellar atmo- 
spheres and molecular absorptions in the spectra of cold 
stars. The Spitzer IRAC camera was partl y designed to 
detect this feature in high re dshift galaxies (|Sawickill2002 : 
Simpson fc Eisenhardt 19991) . The peak is fully sampled 
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by the IRAC photometric bands when it faUs in the 4.5 or 
5.8 /xm channel, i.e. when at least one band lies shortward 
and one longward of the peak. We thus focus our anal- 
ysis on 4.5/xni-peakers and 5.8/xm-peakers (objects with 
SEDs peaking in the 4.5/(ini and 5.8/xm channels). The 1.6 
/zm peak is redshifted to the 4.5 and 5.8 channels for 
sources in the redshift range z = 1.5 — 3. 

We focus this paper on the ELAIS-Sl sub- area cen- 
tered at RA = 00''34"'48" Dec = -43°30'38", where all 
BVR, J, Kg, IRAC and MIPS photometric data are avail- 
able, for a total area of 1 deg^. This area includes 48101 
SWIRE sources, and 29859 SWIRE+ESIS matches. In 
this area, 2848 objects are detected at 24 /im by the MIPS 
camera, while the remaining Spitzer sources are detected 
by IRAC only 

We first apply a 5.8 /zm flux cut at the 3a SWIRE 
depth (25.8 /iJy), in order to favor sources detected in the 
IRAC bands, which sample the restframe near-IR emission 
and hence the stellar mass (as dominated by low mass 
stars). This sub-sample consists of 5546 objects (over one 
square degree). 

3.1. IR-peak galaxies 

The IR-pea k selection is well exemplified in the IRAC 
color space (jLacv et al. l l2004l:IStern et al.ll2005l) . as shown 
in the top panels of Fig. [T] Here red triangles repre- 
sent 5.8/zm-peakers: formally sources showing 5^(3.6) < 
5^(4.5) < 5^(5.8) > S^{8.0). Blue circles (either filled or 
not) are 4.5/xm-peak objects, characterized by 5^(3.6) < 
5,^(4.5) > S',y(5.8). The difference between filled and open 
symbols is defined later in Section [373] We include in the 
plots also objects not detected in the 8.0 /im band, but 
for which the 8.0 /xm 3a upper limit is consistent with the 
definition of the IR-peak. The small dots represent the 
general galaxy population detected in the central ELAIS- 
Sl SWIRE field, light blue having a 4 band detection, and 
orange having a 8.0 /Ltm upper limit only. 

We overplot a set of template tracks as a function of 
redshift: a Seyfert-1 (dotted lines, Mrk 231, Fritz et al., 
20061 ). a Sey fert-2 (dot-dashed, IRAS 19254-7245, Berta 



et al., 2005 ). a starburst galaxy (solid line, M82, Silva 
et al., Il998l en hance d with observed PAHs by Forster- 
Schrei ber et al. 120011 ) . a spiral (long-dashed, M51, Silva 



et al., Il998l ) and an optica lly-blue spiral (short-dashed. 



NGC4490, Silva et al.. ll998h galaxy. The tracks are limited 
to the range z = 1 — 3, for the sake of clarity, apart for 
the starburst one, which is extended down to z ~ 0. 

3.2. Possible aliases 

IR-peak-like colors could in principle be produced not only 
by the 1.6 /xm feature redshifted in the IRAC domain, but 
also by a strong 3.3 fim PAH feature at lower r edshift, as is 



of the IRAC channels, resulting in a blue NIR-IRAC ob- 
served color. 

As far as 5.8/im-peakers are concerned, the require- 
ment that S,y{3.6) < 51,(4.5) automatically avoids low- 
redshift interlopers, because it forces rejection of those 
sources for which the 1.6 /im peak falls at wavelengths 
shorter than the IRAC 3.6 /im channel. 



3.3. Use of near-IR and optical photometry 

On the other hand, in the case of 4.5/im-peakers, only 
one IRAC band samples the SED on the blue side of the 
observed peak and there is no trivial way to avoid aliases 
with low-redshift objects with a bright 3.3 /im PAH. In 
this case the data are of fundamental importance in 
breaking the degeneracy and ruling out interlopers. The 
middle panel in Fig. [T] involving the Ks magnitude shows 
how near-IR data can effectively break this degeneracy 
and help in selecting z = 1 — 3 objects only. The redshift 
tracks of starburst, spiral and blue spiral galaxies suggest 
that objects with {Kg — 3.6) ab < he at z < 1. Therefore 
we adopt this value in order to disentangle low- and high-z 
sources. 

In the two top panels of Fig. [1] (discussed above) , open 
circles represent objects with {Ks — 3.6)ab < 0, while filled 
circles lie above this threshold. In the top left diagram 
( Lacv et al. 20041 ). these low-redshift interlopers seem to 
be separated from the high-rcdshift objects and we empir- 
ically define the transition line: 



(4.5 - 8.0)ab = 1.5 X (3.6 - 5.8)ab + 0.7. 



(1) 



We then use this line to distinguish between low- and high- 
z candidates among those sources that are not detected in 
the near-IR (J, Kg) survey. 

It is also interesting to test whether the optical-IRAC 
color space shows any segregation of low-z versus high-z 
objects (as defined on the basis of the IR-peak technique 
and Ks band photometry), in order to define alternative 
selection criteria to be used when no near-IR data are 
available. The bottom plots in Fig. [T] show a couple of 
diagrams involving th e available optical BVR photometry 
from the ESIS survey ( Berta et al. 20061 ). The majority of 
IR-peakers with {Kg— 3.6) ab < fie in the {Rc— 3.6) ab < 
2 sub-space and below the line 



{Rc - 3.6)ab = 1.3 X {Bj - Rc)ab- 



(2) 



demonstrated by Lonsdale et al. (in prep.) and lBerta et al 



( 20071) . In this case, the real 1.6 /xm peak lies shortward 



In the end, combining the IR-peak criterion and the near- 
IR and optical constraints, our sample contains 149 and 
231 sources peaking at 4.5/tm and 5.8/im respectively, 
over one square degree. Among these, 149 (i.e. ~39%) 
have an optical detection (at least one band), 295 (~78%) 
are detected in either J or Kg, and 109 (^29%) have a 
24/tm counterpart. Finally, ~77% of the IR-peakers have 
no 8.0/im detection. Table [T] summarizes these numbers. 
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(K,-3.e)„ 




Fig. 1. Selection of the IR-peak galaxies. The two top panels show the position of 4 .5Mni-peakers (solid circles) and 
5.8/im-peakers (triangles) in the IRAC color space (jLacv et al.ll2004t IStern et alJuOOST ). The middle plot highlights the 
use of Kg band data to recognize low-redshift interlopers aliasing 4.5/j,ni-peaker colors. The open circles in the other 4 
panels belong to sources with {Kg — 3.6) < [AB mag]. Finally, the two bottom panels present additional information 
based on optical-IRAC colors. Small dots represent the general galaxy population in the SWIRE ELAIS-Sl field. 
Template tracks are overlaid for a Sy-1, Sy-2, starburst (M82), spiral (M51) and blue spiral (NGC4490) galaxies. The 
covered redshift range is reported aside each track. 



Berta S., et al.: IR-peakers mass function 



7 




Fig. 2. Effect of sky-noise on IR-peakers selection. Random fluctuations of IRAC and Kg fluxes can cause objects to 
be scattered in and out of the sample. Solid lines are the actual distribution of IR-peakers. Shaded areas are the result 
of simulations and represent the variation in the number of peakers due to random fluctuations of their colors. The 
three vertical lines refer to the 3, 4, 5 ct thresholds at 5.8/xm in the SWIRE ELAIS-Sl survey. 



Description 


Number 


SWIRE sources 


48101 


SWIRE -1- opt. sources 


29859 


24/xm sources 


2848 


5(5.8) > 25.8 AtJy (3<t) 


5546 


IR-peakers 


380 


IR-peakers 4(t 


326 


4.5/im-peakers 


149 


4.5^m-p. + 24^m 


44 


4.5^m-p. -I- optical 


69 


4.5/im-p. + NIR 


140 


4.5/im-peakers 4(t 


123 


5.8pm-peakers 


231 


5.8/im-p. -I- 24^m 


65 


5.8^m-p. -I- optical 


80 


5.8Atm-p. + NIR 


155 


5.8/im-peakers 4(t 


203 



Table 1. Summary of IR-peaker selection in the central 
square degree of the SWIRE ELAIS-SI field. The descrip- 
tions optical" and NIR" refer to sources with at 
least one optical or JK detection. 



3.4. Effect of sky noise 

Random fluctuations of fluxes within the photometric un- 
certainties in the bands defining the IR-peaker selection 
can cause non-peaker objects to be scattered into the sam- 
ple and actual peakers to fall out of it. 

It is worth noting that, since the 3.6/xm and 4.5/xm 
detection thresholds are much fainter than at 5.8fim (see 
Section[2]), our sources are detected above ^15a in the two 



bluest IRAC channels, and sky-noise affects their fluxes by 
a factor smaller that ^--^5%. 

On the other hand, random fluctuations in the 5.8^m 
and 8.0/im bands can significantly modify the colors of 
these objects and scatter them in and out of the sample. 

In order to quantify this effect, we have performed a 
bootstrap simulation of IR-peaker selection, starting from 
the multi- wavelength SWIRE catalog in the ELAIS-Sl 
field (48101 sources). The procedure applies a random 
fluctuation to the Kj-to-S.O/i fluxes of all SWIRE sources, 
with a dispersion given by the sky-noise associated to each 
source. 

This process was looped 10000 times and a new IR- 
peaker selection was performed at each step. Figure [5] re- 
ports the result of this simulation, for 4.5//m-peakers (left 
panel) and 5.8/im-peakers (right). Solid lines represent the 
actual distribution of IR-peakers, while shaded areas are 
the result of simulations. Vertical dashed lines represent 
the 3, 4, 5 cr detection thresholds in the SWIRE ELAIS-Sl 
5.8/xm survey. 

This analysis shows that, at the 3 a detection level 
random noise fluctuations significantly affect the selection 
of IR-peakers, and tend to diverge. On the other hand at 
the 4 and 5 a flux levels the contamination decreases to 
^ 20% and ^ 15% respectively. 

Therefore the subsequent analysis will be carried out 
only for those IR-peakers detected above the 4 a level (34.4 
/iJy) in the 5.8 ^m band. Thus the sample of IR-peakers 
reduces to 326 objects, 123 of which are 4.5/im-peakers 
and the remaining 203 are 5.8/im-peakers. The electronic 
Table associate to this work reports the main data of the 
final selected sample. 
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4. Photometric redshifts 

The spectroscopic survey in ELAIS-Sl has targeted 
mainly X-ray sources and generally low-redshift (z < 1) 
galaxies. Unfortunately, the few sources detected in the 
IR-peaker redshift range are classified as AGNs, on the 
basis of optical spectroscopy and show power-law IRAC 
colors or bright 8.0/im excesses due to torus warm dust. 
These objects do not belong to the IR-peaker sample. 

As a consequence, we need to rely mostly on photo- 
metric redshift. Nevertheless, recent spectroscopic anal- 
yses of subsamples of IR-peakers in other SWIRE fields 
(Lockman Hole, ELAIS-Nl, ELAIS-N2) have confirmed 
that the adopted selection effecti v ely id entifies galaxies 
at z ~ 1.5 - 3.0. IWeedman erall (|2006l ) performed IRS 



( Houck et al.ll2004 ) mid-IR spectroscopy of IR-peakers de- 
tected at 24/im and found redshifts between z = 1.6 — 2.0 
in 90% of the examined cases. Berta et al. (2007) obtained 
Keck UV-optical (restframe) spectroscopy of IR-peakers, 
confirming the photometric selection between z = \.?> and 
2.5. 

We have d erived photometric red shifts by using the 



Hyper-z code ( Bolzonella et al. 2000l ). We adopt a semi- 



empi rical template library including GRASIL (Silva et al. 
19981 ) models of spiral and elliptical galaxies, M82 
and Arp220 templates (Silva et al.) upgraded with 
observed PAH mid-IR features (Forster Schreiber et al 



I2OOII: ICharmandaris et al.lll999|). and a ULIRG template 
(IRAS 19254-7245. iBerta et ahlboOSl) . 

The Hyper-z performance has been optimized for 
galaxies: we have excluded from the analysis all sources 
with a type-1 AGN spectroscopic classification and all 
those showing a power-law like IRAC SED. In this way 
AGN-dominated objects are avoided. This class is partic- 
ularly delicate to fit with a template-based procedure and 
photometric redshifts are typically mis-int erpret ed in 50% 
of the cases (see for example Berta et al. 20071) . because 
sharp features are missing in their SEDs. Moreover, the 
IR-peak criterion automatically selects the sample against 
power-law AGNs. 

The estimate of photometric redshifts was performed 
including the optical, J, and IRAC (3.6 — 8 /ini) data 
in the computation, but ignoring the MIPS 24/im flux. 
In fact, including the 24/im data turned out to produce a 
higher degree of degeneracy and aliasing. 

The choice of templates in the library and the allowed 
extinction have been tested on galaxy populations in the 
ELAIS-Sl field, taking advantage of the available spec- 
troscopic redshifts. Figure [3] reports the comparison of 
spectroscopic and photometric redshifts, as obtained with 
Hyper-z, The dashed and dotted lines represent ±10% 
and ±20% uncertainties. The difference between the two 
estimates has an r.m.s. of 0.19 and s.i.q.r I of 0.076. 



^ s.i.q.r. = semi inter-quartile range, defined as (53 — ?i)/2, 
where q\ and qz are the first and third quartiles. The first 
quartile is the number below which 25% of the data are found 
and the third quartile is the value above which 25% of the data 
are found. 
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Fig. 3. Comparison of photometric and spectroscopic red- 
shifts for all galaxies with available spectroscopy in the 
ELAIS-Sl field. Objects classified as type-1 AGNs on the 
basis of spectroscopy or having a power-law like IRAC 
SED have been excluded from this analysis. Dashed and 
dotted lines represent ±10% and ±20% uncertainty levels. 
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All IR-peakers 

4.5//m-peakers 
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Fig. 4. Distribution of photometric redshifts of IR- 
peakers in ELAIS-Sl, as derived with the code Hyper-z 
(|Bolzonella et al.ll2000( l. 



The main outliers are sources with few photometric 
data available, especially those with a detection in the 
two more sensitive IRAC channels (3.6/tm and 4.5/tm) and 
one optical band (typically R) only. The photometric code 
calibrated in this way was then used to derive redshifts 
for the general IR-peak population selected as described 
above. 
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The redshift distribution of IR-peakers is shown in Fig. 
m where 4.5^m-peakers are represented by the dotted hne 
and 5.8/xm-peakers by the dashed one. 

5. The estimate of stellar mass 

The estimate of the stehar mass of the galaxies 
in our sample has been obtained with the code 



Sim-Phot-Spec (SPS), dev eloped in Padova (jBerta et al 
2004 IPoggianti et al.ll200lh . 



This code performs mixed stellar population (MSP) 
spectro-photometric synthesis. Several phases in the life of 
a simple stellar population (SSP) are combined together, 
adopting a different star formation rate (SFR) for each 
age. Each SSP phase is meant to represent a formation 
episode of average constant SFR, over a suitable time pe- 
riod Af. The effective number of SSP ages involved in the 
fit depends on the redshift of the given source, in order 
not to exceed the age of the Universe at that redshift. The 
maximum number of SSP phases considered for our sam- 
ple is 6, corresponding to a redshift z = 1, and decreases to 
4 if z = 3. Moreover, during the minimization phase, the 
code checks for SSPs that do not effectively contribute to 
the emitted light, either in the optical and IR, and deletes 
them. In the end the number of SSPs contributing to the 
final fit is typically 3-4. 

Age-selective extinction is applied, assuming that the 
oldest stars have abandoned the dusty medium long ago. 
Keeping in mind that disc populations are on average a f- 
fected by a moderate Ay (< 1 mag, e.g. iKennicutt 1992 ). 
the maximum allowed absorption for stars older than 1 
Gyr is Ay — 0.3 — 1.0 magnitudes. For younger popula- 
tions the color excess gradually increases, but is limited 
to Ay < 5. 

The best fit is found by minimization of the dif- 
ferences between the observed photometric data and the 
synthetic SED, including photometric uncertainties in all 
bands. Photometric red shifts, as computed with Hyper-z 



( Bolzonella et al. 2000L see Sect. HJ are adopted. The 



SPS code takes into full account the uncertainty associ- 
ated with the photometric redshift estimate, while explor- 
ing parameter space. This uncertainty is of the order of 
0.1 — 0.2, depending on how many photometric bands are 
available. When the given source is not detected in one 
(or more) bands, upper limits are used. 

For each galaxy in the sample, the SPS code builds a 
large number of models (~ 10^), exploring the parameter 
space by means of the "Adapt i ve Si r nulated Ann ealing" 
algorithm (ASA, llngbeil l200ll . Il989l : iBerta et all [2004). 
First the parameters are varied with a coarse resolution, 
and then the algorithm focuses on the found minima. In 
order to avoid "freezing" in local minima, re-annealing is 
applied and the code literally "jumps" out of the mini- 
mum in order to explore a different region of the param- 
eter space and find the absolute best fit. In this way, the 
possible fluctuations of colors, within the measured pho- 
tometric errors, are accounted for and are automatically 
propagated to the output stellar mass uncertainty. 



The SPS code is optimized to derive the assembled 
stellar mass in galaxies and the associated uncertainty due 
to degeneracies in star formation history (SFH) space (see 
Berta et al.l 120041 ). The of each model considered dur- 
ing the minimization (i.e. each combination of parameter 
values) is recorded and 1, 2, Scr contours are computed. In 
this way, the resulting uncertainty on the stellar mass esti- 
mate accounts also for minimal- and maximal-mass mod- 
els, derived from the projection of the parameters space 
into "young" and "old" sub-spaces. The resulting best fit 
masses and 3 a mass ranges are reported in an electronic 
Table. 



The 
Padova 



([Fagotto et alj 



adopted SSP hbrary 
evolu ti onary seq uences 
1994allbl: 



IS 

of 



Bressan et al 



based on the 
stel lar models 
1993) 



and 



isochrones ( Bertelli et al.l Il994[ ). and was computed 



by assuming a solar metallicity and a Salpeter initial 
mass function (IMF) between 15 and 120 Mq. The 
SSP spectra were built with the IPickle^ ( 19981 ) spectral 
library , and extended to the near-IR with the iKurucz 
(Il993i ) atmosphere models. Nebular featu res were adde d 
through the ionization code CLOUDY (jPerlandl [l99i ). 
The spectra thus obtained provide a reliable description 
of simple stellar generations up to ~ 5 /xm (restframe). 
Beyond this wavelength, dust emission is no longer 
negligible. 

We assume that the total energy absorbed in the UV- 
optical domain is processed by dust in the thick molecular 
clouds embedding young stars and re-cmittcd in the mid- 
and far-IR (8-1000 /im) in the form of a starburst tem- 
plate. By convolving the template with the 24/xm pass- 
band, the MIPS 24/ini flux expected for the given model 
is computed. Finally, this synthetic flux is compared to 
the observed 24/im data and included in the computa- 
tion. Mid-IR photometry is mainly sensitive to the power 
of the ongoing starburst, as well as to the amount of dust 
obscuring it, therefore it provides a valuable constraint 
on the amount of du st and the streng th of young stellar 
populations (see also lBerta et al.ll2004^ . 

Typically an M82 template is adopted, but a dif- 
ferent choice would not significantly affect the result 



i ng st ellar mass estimate, as demonstrated in lBerta et al 



((20041) . Increasing observational evidence exists that high- 
z IR-peakers detected in the mid-IR resemble the M82 
prototype. Based on Spitzer mid-IR IRS spectroscopy, 
Weedman et al. ( 20061 ) found that z ~ 1.9 IR-peak galax- 
ies are dominated by brig ht PAH features and lack dee 



silicate 10/xm absorption. iRowan-Robinson et al.l ()200 



studied and classified the SEDs of SWIRE sources over 
6.5 deg^ in the ELAIS-Nl fields, finding that M82-like 
starbursts are 3 times more numerous than colder Arp220- 
like objects. Finally, millimeter (250 GHz) observations of 
SWIRE 24^m-selected IR luminous galaxies, performed 
with the MAMBO bolometer on the IRAM/30m telescope 
(Lonsdale et al., in prep.) showed that the 1.2mm/24/im 
flux ratio of these sources resembles that of M82, lower 
that for an Arp220-like population. 
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Fig. 5. Examples of SED fits of IR-peak sources. The top panels belong to 4.5/im peakers, the bottom ones show two 
5.8/ini-peak galaxies. The two sources on the left are not detected at 24/im, and only an upper limit is available, while 
those on the right have an actual 24/xm flux above 250 /iJy (i.e. the 3cr SWIRE limit in ELAIS-Sl). The dashed line 
represents the contribution to the best fit model by young and intermediate-age (<1 Gyr) stars, while the dotted line 
is the light emitted by older stars (age >1 Gyr). The solid line is the total best fit emission up to 5/im (restframe). The 
long-dashed line longward of 5/im (restframe) is the adopted starburst template re-processing the UV-optical light 
absorbed by dust to the IR. 



Figure \5\ shows a few examples of SED fits of IR-peak 
sources. Overplotted on the observed fluxes is the best 
fit solution of the MSP synthesis: the dashed line refers 
to young-intermediate (age < 1 Gyr) stellar populations, 
and the dotted line is the contribution from old (age > 1 
Gyr) stars. The solid line provides the total emitted light 
in the optical and near-IR wavelength range. Longward 
of 5/im (restframe), the long-dashed line is the starburst 



template used to model the IR emission from dust, heated 
by young stars. The open triangle represents the 24/im 
flux predicted by the model, which often is beneath the 
observed data-point. 

Figure [6] reports the distribution of stellar masses as 
a function of redshift. Circles belong to 4.5/im-peakers, 
while triangles represent 5.8/im-peak sources. 
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Fig. 6. Distribution of stellar masses as a function of redshift. Circles and triangles represent 4.5/im- and 5.8/im- 
peakers respectively. Left panel: filled and open symbols belong to sources dominated (best fit) by young or old stellar 
populations. The mass thresholds expected for the maximal and minimal M^/ L ratios in the sample are overplotted 
as dashed (5.8/im-peakers) and dotted (4.5/im-peakers) lines. Histograms are plotted for 4.5/im-peakers (dotted) and 
5.8/im-peakers (dashed). The solid histogram is the total distribution of sources. Right panel: filled and open symbols 
belong to sources with or without a 24/im detection. Histograms have the same meaning, dotted ones representing 
objects detected by MIPS and dashed ones objects not detected. 



In the left panel, open symbols indicate objects whose 
best fit synthetic SEDs are dominated by old (age > 10^ 
yr) stars, while for filled circles the bulk of the luminos- 
ity is provided by younger populations. The dashed and 
dotted histograms belong to 5.8/im- and 4.5/ini-peakers. 
Objects with a resulting lower mass are dominated by 
younger stellar populations, as expected on the basis of 
the dependence of the Mi,/L ratio on age. The effect of 
Malmquist bias is evident not only from the trend seen in 
the Mi, vs. z plot, but also from the distribution of sources 
in the mass space: 4.5/im-peak galaxies have masses lower 
than 5.8/im-peakers, on average. 

Typical error bars are shown: uncertainties on the stel- 
lar mass can be as high as 0.3 dex for sources that lack two 
{B and V) or all optical bands and — in the worst cases 
— also near-IR photometric data. This situation is repre- 
sented by the big error bar. When good multi-wavelength 
coverage is available, the uncertainty on the stellar mass 
reduces significantly (small error bar). 

The mass thresholds expected at the 4 a flux limit 
of 5'(5.8/im) = 34.4 fiJy, as derived from the maximal 
and minimal Mj,/L ratios in the SED fitting results, are 
overlaid on the data. Dotted lines belong to 4.5/im-peakers 
and dashed fines to 5.8/im-peakers. See Sect. 16.21 for a 
discussion on the completeness of the sample. 

The right panel of Fig. [6] shows the same plot, but 
coded on the basis of 24/im detections. Filled and open 
symbols refer to sources detected or not detected in the 
MIPS 24/im channel, respectively, at the 3cr level. The 
dashed and dotted histograms have the same meaning. 



5.1. Effect of 24 /im data on the mass 

The availability of mid-IR observations provides a valuable 
tool for constraining the amount of dust and the luminos- 
ity of young stellar populations in the synthesized models, 
thus reducing the unce rtainty on the derived stellar mass 



(e.g. iBerta et al.ll2004l) 



In order to understand the effect of the 24/im flux on 
the mass, i.e. how it effectively influences the star forma- 
tion history of the best fit solution, we have performed a 
fitting run ignoring the available 24/im photometry. 

The results are shown in Fig. [T] The top-left panel 
shows the distribution of stellar masses as a function of 
redshift (symbols are the same as in Fig. [6]), while the 
top-right plot compares the new results (solid line, ob- 
tained excluding the 24 /im fiux from the minimization) 
to the standard fit mass distribution (dashed histogram). 
The two bottom panels present the same comparison, for 
sources dominated by young and old stellar populations 
respectively (more that 50% of the total mass being as- 
sembled in stars younger or older than 1 Gyr). 

If no mid-IR detection is available, the amount of dust 
and young stars are free to vary with no control, and 
very extinguished [Ay > 4) young populations, not vis- 
ible in the optical domain, may exist. As a consequence, 
the amount of young stellar mass is found to be larger, the 
total stellar mass inferred from SED fitting smaller on av- 
erage, and the mass spread slightly wider. The availability 
of MIPS data provides a valuable constraint on the recent 
history of star formation of galaxies and helps to avoid a 
significant number of SED fitting solutions that could not 
be a priori ruled out. 
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Fig. 7. Results of SED fitting, ignoring the 24/im fluxes (solid histogram), and comparison to the fit including 24^m 
(dashed, see Fig. [S]). In the top left panel, filled and open symbols represent sources dominated (best fit) by young 
or old stellar populations, respectively. The mass distributions for these two sub-classes are shown in the two bottom 
panels. 



5.2. Effect of near-IR (JK) data on the fit 

The near-IR J and Ks data turned out to be very useful 
at the selection stage, in particular to exclude low-redshift 
interlopers with IRAC colors similar to 4.5/xm-peakers. 

It is worthwhile exploring how observations at these 
wavelengths (1-2 /im) constrain the estimate of stellar 
masses, by means of SED fitting. We have therefore per- 
formed another fitting run, without taking into account 
the J, Ks photometry. Figure [8] shows a couple of SEDs 
and the comparison to masses from the standard fit. 

The top panels include the SEDs of a 4.5/im-peaker 
(left) and a 5.8/Ltm-peak galaxy (right). The dashed fine 
represents the fit obtained with the standard technique, 
i.e. accounting for all the available photometric data, while 
the dotted lines are the best fit models obtained by exclud- 
ing the near-IR (J and Ks) fluxes from the minimization. 

In some cases, the J, Ks data turn out to be use- 
ful to constrain the age of the dominant stellar popula- 
tions hosted by IR-peakers, since they sample the depth 
of the D4000 break, when combined with optical data. 
Nevertheless, the J, Ks constraint turns out to only be 
effective, in this respect, for a small fraction of cases 
(< 15%). The top panels in Fig. [5] show two examples 
among those for which these data are most useful. 



It is worth noting that the advantage of having near-IR 
data is larger for 4.5/im-peakers than for objects that peak 
at longer wavelength (in the 5.8 /im channel, i.e. at higher 
redshift). For the latter, the SED slope blueward of the 
1.6/im peak is, in fact, already deflned by the 3.6—4.5 — 5.8 
fim colors well enough to reasonably constrain the D4000 
break, and therefore the best fit solutions obtained with 
or without J, Ks do not differ too much. 

Several authors fe.g. lBerta et al ]L2004') had shown that 
the introduction of IRAC data in SED fitting provides 
tighter constraints on the stellar mass, reducing its un- 
certainty b^;_a^Jactor as high as 5, for z > 2 sources. In 
addition, Wuyts et alJ ()2006,) infer that JHK band data 
can reduce the uncertainty on the stellar mass of blue 
galaxies with {U — t^)i-est < Ij while IRAC photometry is 
really effective only for redder sources, having a _R— IRAC 
color significantly redder than 1; this is the case for our 
IR-peakers. 

Overall, the distributions of stellar masses, as obtained 
with and without the near-IR J and Ks data (Fig. [5]), are 
not significantly different, confirming the result that these 
data do not play a fundamental role in deriving the stellar 
mass of IR-peakers. On the other hand, as we have already 
pointed out, near-IR magnitudes are critical to avoid low- 
z aliases. We also expect these data to be more effective in 
constraining the age and M^/L ratio of galaxies in other 
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redshift bins, e.g. when the blue side of the 1.6/im feature 
is not fully sampled by IRAC, z < 1.5. 



5.3. Choosing the IMF 

The choice of an initial mass function (IMF) different from 
Salpeter's can significantly affect the estimate of the stel- 
lar mass. 

The lSalpeteij (Il955[ ) IMF has a constant a = 1.35 slope 
throughout the whole mass range considered, but in fact 
it was never measured down to 0.15 M0 by Salpeter. 



More recent determinations of the IMF showed that a 
flatter slope is needed at low masses (M < IM^, in order 
to reproduce the observed Galactic data (IMiller & Scalo, 



1979 : Kennicuttl[l98l iKroupa et aTlll993inKimrF3l200ll 



Chabrien .2003al . e.g.) Consequently, assuming Salpeter's 
IMF results in including too many low-mass stars (which 
dominate the stellar mass budget) in the galaxy modeling. 
Introducing a drop-off at M < 1 Mq leads to lower values 
of the stellar mass of the analyzed galaxies. 

At the bright end, lElmegree 3 (I2OO6I) pointed out that 
above 1 Mp, the IMF slope is not steeper than Salpeter's; 
Miller fc Scalol (|l979l ) and the other IMFs with a drop at 
the highest masses were based on galactic disk measure- 
ments which cannot be safely used to trace the high-mass 
end of the IMF, because of the complicated SFH of the 
Galaxy. 



We study here the effect of a different choi ce in the 
IMF, b y performing a new fitting run with the Chabrier 
(|2003ah IMF instead of the classic [Salcete] (|l955l ) one. 



The IChabrieij (j2003a[ ) IMF is described by a power- 
law for M > 1 M c^ an d a logno r mal fo rm below this limit. 
Chabrieij (|2003ah and iKroupal (|200lh are very similar to 



each other, but here we prefer to adopt the former be- 
cause it is physically motivated and provides a better fit to 
counts of low-mass stars and b rown dwarfs in the Galactic 
disc (see for examp le Ch abrier l200l|; |200i [2003^ and also 
Bruzual & Chariot llooi). 

The results of this analysis are reported in Fig. [9l 
where we show the shape of different IMF's (top left 
panel), and the direct comparison between stellar masses 
as obtained with the Salpeter and Chabrier IMFs (top 
right panel). In the bottom panels, the comparison of 
the stellar mass distribution for the two cases is shown. 
Solid histograms represent stellar masses derived with the 
Chabrier IMF and dashed ones with Salpeter's. 



The diagrams show that the IChabrieil (l2003aD IMF 
leads to systematically lower masses than the ISalpeter 
(jl955l ). as expected. The difference between the two turns 
out to be '^0.3 dex, represented by the dashed line in the 
top right panel of Fig. [9l The solid line sets the 1:1 ratio. 

The stellar mass distribution is roughly rigidly shifted 
to lower masses (bottom panels). As far as the splitting 
into sources dominated by young-intermediate (age < 10^ 
yr) and old stars (> 10^ yr), a few objects migrate from 
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the former to the latter sub-samples, but the overall dis- 
tributions are maintained. 



Despite t he fac t that the choice of a IChabrieij (|2003a[ ) 



or 



Kroupal (|200lh IMF restricts the inclusion of too 
many low-mass stars, the majority of literature studies 
on the stellar m ass function of galaxies ar e based on the 
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1998f ). For this reason, we will derive the stellar mass 
function of IR-peakers from the Salpeter-based SED 
fitting, keeping in mind that a different choice (e.g. 
Chabrien l2003ar ) would produce a shift of the mass 



function to lower masses (e.g. by about 0.3 dex). 
5.4. Models with TP-AGB phase 

Recently, iMaraston et al. I (|2006h have performed stellar 
population synthesis of z = 1.4 — 2.7 galaxie s from the 
GO ODS Spitzer surv ey ( Dickinson et al.ll2003a[ l. adopting 
the MarastonI ( 2005h library and studying the differences 



in the derived par ani eters, with respect to the Padova (e.g. 
Fagotto et al. 1994a) models. 



The iMaraston (■2005i') library is based on the Frascati 
tracks (e.g. Cassisi et al.lll997h . The main difference be- 



tween the two is that the Padova tracks include a certain 
amount of convective overshooting on the main sequence 
(MS), whereas the Frascati tracks do not; moreover the 
temperature distribution of the red giant branch (RGB) 
phase is shifted to cooler temperatures in the Padova 
tracks for Z > Zq. As a consequence, the MS lifetime 
is longer for the Padova tracks and the RGB phase is de- 
layed, with respect to the Frascati models^ 

In any case, the key difference of the MarastonI ( 2005 ) 
approach is the way the thermally-pulsing asymptotic gi- 
ant branch (TP-AGB) phase is included in the evolution of 
stellar populations, i.e. by means of a semi-empirical fuel 
consumptions table, in contrast to "a posteriori" recipes 
used in isochrone synthesis. TP-AGB stars dominate the 
near-IR (e.g. K band) luminosity for SSP ages between 
0.3 and 2 Gyr (for a Salpeter IMF and solar metallicity), 
while main-sequence stars dominate in the optical (e.g. V 
band) . 

At the transition between the early-AGB branch and 
the TP-AGB, i.e. when the TP-AGB phase begins, the 
near-IR luminosity of stars significantly increases, and 
therefore the near-IR mass-to-light ratio (M^/L^f) of the 
SSP drops suddenly by a factor 3 — 5. In models without 
this phase, M*/L monotonically increases, independent of 
wavelength. At ages between 0.8 and 1 Gyr, the models 
based on the Padova tracks have M^^/La' larger than the 
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Maraston ones, as a consequence of the different 



treatment of the TP-AGB phase. At older ages, M^^/Lk 
is smaller for the Padova 1994 case, because of the cooler 
R GB phase. 



Maraston et al 



(|2006l ) model the observed SEDs of 
high-z galaxies, adopting different star formation histories 
(SFH): an instantaneous burst, a exponentially declining 
star formation (SF) , a prolonged b urst, and a const ant 
SF. As a result, they find that the iMarastonI (|2005[ ) li- 
brary leads to systematically younger best fit solutions 
for these systems, and hence to lower stellar masses, with 
respect to the Padova 1994 library. The IRAC Spitzer data 
tur n out to be very useful in constraining the mass, and 
the iMarastonI (j2005l ) models provide a better fit to the 
observed photometry. 

With this in mind, we have performed a si milar anal- 
ysis o n our sample of IR-peakers, adapting the Marastoiil 
Hooi) SSP library to the SPS code and running it in the 
same exact way as before. Figure[Tn]compares the SEDs of 
the Padova and Maraston SSPs, in the age range in which 
the TP-AGB phase is active. 

The results of this further analysis are reported in Fig. 
[TTl where the stellar masses derived with the Maraston 
library (solid histogram) are compared to those obtained 
earlier with our "standard" fit (dashed) . The top-left panel 
reports the M^/L^ q ratio for the SSPs in the two libraries 
and the top-right plot shows the direct comparison of stel- 
lar masses as derived in the two cases. As far as the low- 
mass end of the distribution is concerned, stellar masses 



in the Padova-94 case. Conversely, the Padova-94 fit pro- 
vides slightly lower masses at the bright end. 

The free-form MSP technique adopted for fitting the 
I R-peakers produces a d ifferent result, compared to that 
of Maraston et al. ( 2006f ) . The new library of models pro- 
duces best fit solutions characterized by completely differ- 
ent SFH's, with respect to the previous ones, based on the 
Padova tracks. The relative contribution of young and old 
stars in the best fits changes significantly, with the frac- 
tion of objects dominated by old populations increasing 
in the Maraston case. The two bottom plots show the M^, 
distribution for sources dominated (in mass) by old and 
young stars: not only the balance between old and young 
SSPs changes, but also the distributions peak at lower 
masses (for old-dominated sources) and higher masses (for 
young-dominated ones), than in the Padova-94 case. On 
the other hand, the overall stellar mass distribution (top 
right panel) is not significantly changed, apart for a shift 
of M < 10^^ M0 galaxies to higher masses. 

5.5. Possible AGN contributions 

First of all, the sample of IR-peakers has been checked 
against bright X-ray sources, th anks to the ayailabl e 
XMM-Newton survey in the area (jPuccetti et al.ll2006l ). 



Ten sources turn out to host an AGN X-ray component, 
seven 4.5(Um-peakers and three 5.8/im-peakers. The for- 
mer are characterized by a bright 8.0/zm excess, likely due 
to the presence o f toru s warm dust in the mid-IR SED (see 
also Berta et al. 20071 Lonsdale et al. 2007, for an exten- 
sive description of SED shapes). The three 5.8/^m-peakers 
emitting in the X-rays show a smooth stellar peak, diluted 
by the AGN dust. These sources have been excluded form 
further analyses. 

The remaining sources show a sharp IR-peak, and, in 
performing the SED analysis, we have assumed that only 
stellar emission contributes to the observed SEDs, while 
no AGN component is present. In fact, the warm dust in 
the AGN torus would emit at IRAC-MIPS wavelengths, 
producing a power-law SED or diluting the 1.6/im (rest- 
frame) stellar peak. Hence a sharp stellar peak shifted to 
the 4.5 or 5.8 /im bands should identify sources whose 
near-IR emission is dominated by stars. Nevertheless, the 
presence of a possible AGN component cannot be fully 
excluded. 

have analyzed UV-optical (rest- 



iBerta et al 



(|2007D 

frame) spectra of IR-peakers observed with the Keck- 
I 10m telescope, selected in order to show no evidence 
of AGN c omponen ts on IRAC c olor-color diagrams (e.g. 
Stern et al . 2005; Lacv et al.ll2004) . Because of instrument 
limitations, the sample was limited to the brightest IR- 
peakers in the SWIRE survey, with r' <^ 24.5 (Vega). 

These authors find evidence for AGN emission lines in 
62% of the detected IR-peakers, two thirds belonging to 
the type-1 and one third to the type-2 classes. The spectro- 
scopic redshift of these sources lies between z = 1.3 — 2.0, 



based on the Maraston library turn out to be higher than including 5.8/xm-peakers. In order to explain the pres- 
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ence of the 1.6/im peak in the 4.5 and 5.8^m channels, 
a multi-component (stars and AGN) fit to the SEDs of 
these sources was performed. The AGN component par- 
tially dilutes the 1.6/im peak, and — more importantly — 
modifies its shape, resulting in an apparent shift to longer 
wavelengths. 

At optical restframe wavelengths, the AGN emission 
is overwhelmed by stars in the host galaxy, while it re- 
emerges in the IRAC domain, especially in the two longest 
wavelength channels (5.8 and 8.0 iim): torus dust is not 
negligible for these sources, and thus might affect the 
1.6/im-peak selection. All the sources hosting an AGN 
were detected at 24 /im above 250 /iJy, have moderate 
24/xm excess ([3.6 — 24] = 1.5 — 2.5 in AB units) and show 
a wide range of optical-IRAC colors (Re — 3.6)^^ = 2 — 4. 
AGN features were identified only in sources above r' <^ 
23.8, corresponding t o Vj ~ 23.6 (AB) for a gray source 



(|Fukugita et al.lll996l ) 



If an AGN component were present, the estimate of 
the stellar mass would be affected in two different ways, 
operating in the same direction, i.e. decreasing the ac- 
tual mass, with respect to that measured by ignoring the 
AGN. Firstly, the torus warm dust emission would re- 
sult in a lower redshift than expected for IR-peakers, and 
the luminosity (hence the mass) would be consequently 
lower. Moreover the torus would contribute to the ob- 



served IRAC fluxes, therefore the light in stars (hence the 
mass) would be even smaller. 

We can not a priori rule out the presence of an AGN 
component; nevertheless, it is worth noting that only 48 
sources (i.e. '^12.5% of the IR-peaker sample in ELAIS- 
Sl) are brighter than Vj — 23.6 [AB]. Among these, only 
23 are also detected at 24/im, corresponding to '^6% of 
the sample, and finally only ~3.5% (13 objects) show a 
moderate (3.6 — 24) color. 

We will not discuss this problem further, since the frac- 
tion of sources that might be affected is very small. 



6. The mass function 

The results obtained from spectro-photometric synthesis 
have been exploited to derive the contribution of z ~ 1 — 3 
IR-peakers to the galaxy global stellar mass density, p*. 

The derived stellar masses are used to build the 
stellar mass function of IR-peakers in t wo ways: wit h 



the V g technique and applying the STY (jSandage et al 
Il979f ) method. The latter uses a parametric function, 
^ (M, param) and a Monte Carlo Markov Chain algorithm 
adopting a Bayesian formalism. 
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on the technique adopted to correct incompleteness. 



6.1. The accessible volume 

In order to compute the comoving number density of IR- 
peakers at the given reds hift, we have a. dopted the well- 
known l/Va method (e.g. ISchmidtl 19681 ). The accessible 
volume, Va, in which each galaxy could be detected in the 
survey is computed by taking into account the selection 
criteria: 

— the flux cut at ^(S.S/im) = 34.4 /iJy, 

— the IR-peak color selection. 



The f ormer provides the classical Vmax estimate ([Schmidt 
1968() , related to the maximum redshift Zmax at which the 



given galaxy would be observable in our survey. The z^ax 
is computed by adopting the best fit SED model (see Sect. 
[5]) , redshifted and k-correctcd until cosmo logical dimm ing 
fades it below the adopted fiux limit (e.g. lHog3ll999f ): 



(3) 



where (11 is the luminosity distance and K(iy, z) is the 
fc-correction through the filter T{v): 



K{v, z) = (1 + z) 



r^lL[v{l + z)]T{v)dv 
rL{^)T{v)di. 



(4) 



The IR-peak condition defines a spherical shell, de- 
limited by redshifts z^^^ and z^^^, where the IR-peak 
selection is valid, i.e. where the 1.6/^m restframe peak is 



detected in the IRAC 4.5 or 5.8/zm channel. The effective 
accessible volume for each source in our color-selected sur- 
vey is thus given by: 



Ait ,p<!<^k az 



(5) 



where VL is the surveyed sky area, and the volume element 
dV depends on the adopted cosmology. 

It is worth noting that the {Kg — 3.6) > cut, applied 
in order to avoid low-redshift interlopers, does not affect 
the accessible volume estimate, because the cutoff redshift 
is always smaller than z^ 



peak 
mm ■ 



6.2. Completeness 

The IR-peaker sample was first selected by applying a 
flux cut at 5.8/im. Despite the fact that we are directly 
probing the restframe near-IR emission of these galaxies, 
which is primarily powered by low-mass stars dominating 
their stellar mass budget, it is not possible to define a 
sharp mass limit encompassing the whole sample. 

In fact, the mass-to-light ratios (M-^jV) of these galax- 
ies spans a relatively wide range, between ~0.03 and ^0.5 
at 3.6/Ltm (restframe). At a given redshift and flux, these 
translate into very different mass values, as shown by the 
minimal and maximal galaxy tracks shown in Fig. [5] (left 
panel). As an example, it would thus be more correct 
to say that at 2; = 2 the sample is limited to masses 
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M > 7 X 10^" or M > 3.5 x 10" Mq for fow or high 
M^/L objects, respectively. 

This effect has been thoroughly ex amine d 



by 

Dickinson et~all (l2003b[) . iFontana et all (|2006l ) and 



Fontana et al. " |2004[ l. The latter perform a very detailed 



analysis and provide a valuable recipe to partially correct 
the mass incompleteness of flux-limited samples. 

Figure [12] reports the distribution of stellar masses as 
a function of the observed 5.8/im flux for the galaxies 
in our sample, split into different rcdshift sub-bins. The 
left panels refer to 4.5/im-peakers, while on the right side 
5.8/im-peak galaxies are shown. The vertical dashed line 
represents the adopted 5.8/xm flux 4 a cut. The diagonal 
dashed and dotted lines are the tracks described by the 
sources with minimal and maximal M-^/ L ratios in the 
sample (see also Fig. |6]), at the central redshift of each 
sub-bin. 

Since the majority of 4.5/im- and 5.8/im-peak sources 
are in the 1.50 < z < 1.75 and 2.25 < z < 2.50 sub-bins 
(see the redshift distribution in Fig. [6]), we will use these 
as reference to derive the threshold completeness mass for 
the two different classes. 

The horizontal lines in the two panels focused on these 
redshift sub-bins represent the stellar mass levels above 
which the samples are definitely complete, because these 
thresholds are set by the maximal M-^jL observed values. 
These masses turn out to be 2 x 10^^ and 3.5 x 10"^-^ 
for the two populations. Below these values, the sample 
becomes incr easingly incomplete. 

Following iFontana et al.l (|2004[) . it is possible to com- 
pute the fraction of sources lost at each given mass, by 
exploiting the observed distribution of M^/L ratios and 
fluxes. We defer the reader to their paper for further de- 
tails on this technique. We keep those mass bins for which 
the fraction of missing sources turns out to be <^ 0.5. In 
this way it was possible to extend our analysis down to 
1.25 x 10" and 1.6 x 10" Mq for ^.h^iTOr and 5.8//m- 
peakers respectively (see Tab. [2]). Note also that below 
these limits the amplitude of random flux fluctuations 
due to sky-noise is large and the contamination of the 
IR-peaker sample by other classes of sources can not be 
controlled (see Sect. 



6.3. The observed galaxy comoving number density 

Within the Va formalism, the comoving number of galaxies 
per unit volume, in each redshift bin and in the mass bin 
AM, is obtained by: 



$(M)AM = ^i-AM, 



(6) 



where the sum is made over all galaxies in the given mass 
bin. 

Figure [T3] shows the distribution of the comoving 
number of galaxies as a function of stellar mass (i.e. 
the "observed" stellar mass function) for 4.5/L(m-peakers 
(dashed line, filled circles) and 5.8/zm-peakers (solid line. 



triangles). Table [5] reports the data. The stellar masses 
come from the spectro-photometric fit obtained with the 
Padova-94 library and the Salpeter IMF, and accounting 
for all the available photometric data and upper limits. 
Filled symbols represent the mass function corrected for 
incompleteness, while open symbols show the data as ob- 
tained before applying any correction. 

The number density of 4.5/im-peak population turns 
out to be significantly lower than that provided by 5.8/L(m- 
peakers. It is worth noting that the selection based on a 
5.8/im flux cut is optimized for the detection of the 1.6/im 
stellar peak in IRAC channel 3, while in order to have 
a comparable selection for 4.5/im-peakers we would have 
needed to apply a 4.5/im flux cut. The 5.8/im flux cut 
corresponds to a restframe H band selection for 5.8/im- 
peakers and to a band selection for 4.5/xm-peakers. 
Assuming that the typical S{H) / S{K) flux ratio of a 
galaxy (e.g. a IR-peaker) is ^ 1.4, then the S'(5.8/xm) = 
34.4 /iJy cut corresponds to S'(4.5/im) = 48.2 /xJy for a 
4.5/im-peak galaxy. As a confirmation of this effect. Figure 
1141 shows the distribution of 4.5/im and 5.8/im fluxes for 
our sources. 

However, the 5.8/im-peakers sample includes also all 
those galaxies not detected at 8.0/im, but still consistent 
with the 5.8/im-peak selection, when using the S.O/xm up- 
per limit. Unfortunately, since the 5.8/im band is the least 
sensitive among the 3.6, 4.5 and 5.8 /im SWIRE/IRAC 
channels, it is not possible to perform a comparable 4.5/im 
selection. In fact, in order to include similar objects to 
the 4.5/im-peakers sample, all the z = 1 — 2 galaxies with 
S'(4.5/im) ~ 34.4 /iJy, not detected at 5.8um, but still 
consistent with the 4.5um-peak selection, should be taken 
into account. The resulting sources are detected only in 
two IRAC channels (3.6/im and 4.5/im) and are barely 
detected in J or Ks. Hence aliasing by low-z interlopers 
or contamination by power-law objects would be hard to 
control. 

We conclude that our completeness correction is not 
able to control this selection deficit in the 4.5/im-peakers 
sub-sample. Therefore we will limit the parametric deriva- 
tion of the mass function to the 5.8/im peakers only. 

The error bars on the mass function plotted in Fig. [12] 
account only for Poisson statistics. Nevertheless, the un- 
certainties on the stellar mass estimate, caused by degen- 
eracies in the SFH space and errors on photometric red- 
shifts, introduce a non negligible contribution to the un- 
certainty on the observed stellar mass function and must 
be taken into account. 

Due to these uncertainties, an object can move from 
the mass bin it formally belongs to and actually contribute 
to the stellar mass function in another mass regime. The 
probability for each galaxy to contribute to the comov- 
ing number density in other mass bins is given by the 
distribution obtained during SED fitting with the ASA 
algorithm. 

For sake of clarity we do not plot these additional error 
bars in Fig. [131 but they will be fully taken into account 
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Fig. 13. Mass function of 4.5/im-peakers (dashed line, circles) and 5.8/im-peakers (solid line, triangles) in the ELAIS- 
Sl area, as obtained with the standard spectro-photometric fit and the Padova-94 library (see Sect. [5]). Filled symbols 
belong to the completeness-corrected mass function, while open ones represent the data before any correction was 
applied. Error bars account only for Poisson noise statistics. 
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Table 2. Observed mass function of SWIRE IR-peakers, as obtained after correction for incompleteness effects, for a 
Salpeter IMF and Padova-94 stellar tracks. We report the measure values of the galaxy comoving number density and 
the fraction of missing sources in each mass bin (see the discussion on completeness in Sect. Wl2\i . For 5.8/im-peakers 
the mass function obtained after evolutionary correction is also included (column 4). 



when fitting the observed data with a parametric mass 
function, by using a Bayesian approach (see Sect. 16. 6p . 

6.4. Results with the TP-AGB enhanced phase 

The IR-peakers mass function obtained with t he Padova- 
94 lib rary is compared to the results from the iMarastonl 
approach in Fig. [151 Data are corrected for incom- 
pleteness, in the same way described for the Padova-94 
case. 

As far as 4.5/im-peakers are concerned, the two li- 
braries lead to very different mass functions, mainly be- 
cause low -mass obj e cts m igrate to higher mass values, 
when the Maraston SSPs are used (see also Fig. 

[TTjand Sect. 15. 4p . Nevertheless, this mass function is still 



affected by a strong incompleteness effect, which cannot 
be recovered. 

In the case of 5.8/im-peakers, the difference between 
the two results is less dramatic and the two libraries lead 
to very similar mass functions. 

Note, however, that the comparison of our results to 
previous estimate s of the stellar mass function at high 
redshift (e.g.iFontana et al.l l2006L [iooi lOrorv et al.ll2005: 
Franceschini et al. 20061 among others), will make use of 
the results obtained exploiting the Padova-94 library. 

6.5. Evolutionary effects 

The mass density of galaxies increases very rapidly be- 
tween redshift 3 and 2, by a factor of ^ 4, according to 
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Fig. 16. Left panel: mass function of 5.8/zm-peakers split into three redshift sub-bins (same symbols as in Fig. [T5)) . 
Right panel: evolutionary trend of the 5.8/Ltm-peakers stellar mass density as a function of redshift for three different 
mass bins. Solid lines connect the observed data, dotted lines represent a power law fit to the evolution. 
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Fig. 14. Distribution of 4.5/im and 5.8/im fluxes for IR- 
peakers. Circles represent 4.5/i-peakers, triangles belong 
to 5.8/im-peakers. The vertical dashed line corresponds to 
the S'(5.8/im) = 34.4 /iJy flux cut. The horizontal dashed 
line sets the 5(4.5/im) = 48.2 /iJy limit, obtained by as- 
suming a restframe S{H)/ S{K) = 1.4 ratio for typical 
IR-peak galaxies. 



Fontana et al. see also Sect. 17. 2| ). It is thus worth- 

while studying the evolutionary details of the galaxy stel- 
lar mass assembly over this redshift range. The left panel 
in Fig. [16] reports the stellar mass function of 5.8/Ltm- 



s 



s 
■a 



10- 



10-s 



10- 



4.5/Lim-peakers :: 5.8ju.m-peakers : 




-•- Padova94 
■o- Maraston05 

j_lJ 




Padova94 
MarastonOS 



10" 5x10" 10" 10'= 

Mass [h.;= Mg] Mass [h;= M^,] 



Fig. 15. Comparison of the mass function as obtained 
with the Pado va-94 (filled symbols, solid lines) and 
Maraston (2005^ open symbols, dotted lines) libraries, af- 
ter correction of incompleteness. 



peakers split for the first time into three redshift sub-bins: 
2.0 < z < 2.25, 2.25 < z < 2.5 and 2.5 < z < 3.0. Note 
that the integration boundaries in Eq. [5] are now set by 
the used redshift sub-bins. 

The number of massive galaxies decreases in the higher 
redshift sub-bins, as the normalization of the MF becomes 
smaller by a factor 2-4, depending on mass. Despite the 
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Mass range 


a 


h 


[10" Mq] 






1.6-3.2 


-0.64 ± 0.06 


-2.27 ± 0.19 


3.2-5.0 


-0.39 ± 0.04 


-3.40 ± 0.09 


5.0-7.0 


-0.10 ± 0.11 


-4.73 ± 0.28 



Table 3. Evolutionary coefficients describing the depen- 
dence of the stellar mass density on redshift and mass, in 
the form of a (1 -f z) power law (see text for more details). 
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Fig. 17. Effect of evolutionary correction on the stellar 
mass function of IR-peakers. The open symbols and the 
dotted line belong to the mass function derived without 
evolutionary correction, while the filled symbols and the 
solid line account for it. 

large error bars, a flattening of the stellar mass function 
seems to be detected, as redshift increases. The right panel 
in Fig. [in] shows the evolution of the stellar mass density 
as a function of redshift and for three different mass bins, 
highlighting that the highest-mass tail (M > 4 x 10" M©) 
of the mass function evolves less rapidly than lower-mass 
bins. 

For each mass range considered, we have reproduced 
the evolutionary trends with a power law: 



dN 
dM 



(z) oc b (M) X (1 -f zf'^^'^^ 



(7) 



The resulting evolutionary parameters are reported in 
Tab. [21 Dotted lines in the right panel of Fig. [TBI represent 
the power law fit. 

In any case, it is worth to point out that above z — 
2.5 too few galaxies are included in the 5.8/im sample, 
and the significance of the derived mass function is rather 
poor. Only by exploiting the whole SWIRE area (49 deg^) 
will it be possible to build a catalog of massive galaxies 
at z = 2 — 3 large enough to study the evolution of the 
number density at these epochs. 



In order to compare our results to l iterature data (e.g. 
Fontana et 31112009: lOrorv et al.|[200l iGwvn fc Hartwick 



2005[ ). the whole z = 2 — 3 bin will be considered in the 
following analyses. Computing the stellar mass density on 
this wider redshift range allows one to exploit a higher 
S/N ratio, but on the other hand dilutes the information 
on the evolution of the mass function across redshift. 

As far as a broad redshift range is considered, it is 
very important to correct the derived stellar mass density 
against evolutionary effects. 

Since this is a flux-limited sample, fainter (less mas- 
sive) objects probe smaller volumes than more massive 
objects. Therefore the highest mass tail of the mass func- 
tion is averaged over the full z range considered, but the 
low-mass bins are contributed by galaxies populating only 
the low-z part of the whole z = 2 — 3 redshift range. Hence 
at the faint-end the effect of negative evolution is not bal- 
anced and tends to steepen the mass function. 

If on one hand the volume effect is corrected by the Va 
and completeness analyses, on the other hand the intrinsic 
evolution is not accounted for in this way. Moreover this 
evolution strongly depends on mass. 

We therefore computed the necessary correction fac- 
tor for each galaxy, by comparing the average stellar mass 
density over the z = 2 — 3 range to that predicted at 
the given redshift by our power-law fit, as a function of 
mass. This additional coefficient is then applied to Eq. [6| 
when computing the stellar mass function. Figure [171 com- 
pares the stellar mass function obtained without and with 
the evolutionary correction. The low-mass end flattens, as 
expected. The fourth column of Tab. [H reports the mass 
function of 5.8/i-peakers corrected for evolution. 

6.6. Parametric stellar mass function 

An independent characterization of the stellar mass func- 
tion of massive galaxies at high redshift is the pa ramet- 
ric approach by Sandage, Tammann, Yahil (Il979l . STY). 
We focus the analysis on 5.8/xm-peakers only, because the 
4.5/xm-peaker sample turned out to be affected by incom- 
pleteness effects that we could not control (see Sect. 16. 3|) . 

The results of this analysis will then also be used to 
estimate the contribution of this population to the global 
stellar mass density of galaxies. 

As commonly found in the literature, we choose to re- 
produce the observed d ata w ith the usual parametric de- 
scription bv lSchechter ( 197d ): 



$ (M) = $* ( — 



(8) 



where — for the sake of clarity — M here denotes the 
galaxy stellar mass, previously called M*. 

The three free parameters $*, a and M* represent the 
normalization, the slope in the low mass regime, and the 
transition mass between a power-law and the exponential 
drop-off or the e-folding mass of the latter. 

Since we are sampling the very high mass tail of the 
galaxy stellar mass function, SWIRE IR-peakers are well 
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suited to constrain the M* parameter, while we expect 
to probe a only in a poor way. We will discuss this issue 
later. 

While performing this analysis, it is very important to 
account for the main sources of uncertainty on the galaxy 
number density, as described above (see Sect. 16. 3p . The 
first source of uncertainty is represented by Poisson noise, 
since the sample consists only of 203 5.8/xm-peakers. The 
second is that the stellar mass estimate for each individ- 
ual galaxy has non-negligible uncertainty, which must be 
accurately included in the parametric fit. 

We used a Markov-Chain Monte Carlo (MCMC) sam- 
pling of the parameter space, to explore the posterior 
probability function of the model, with parameters com- 
prising both the three Schechter parameters and also ad- 
ditional hyper-parameters to represent each of the indi- 
vidual galaxy masses. 

According to the Bayes' theorem: 



v{e\d) 



r{d\e)v{e) 
v{d) ■ 



(9) 



i.e. the conditional probability of the Schechter parame- 
ters, given the data, is equal to the conditional probability 
of the data, given the parameters, times the prior proba- 
bility of the parameters, divided by a normalization factor. 
The theorem can be paraphrased as: 



posterior = 



likelihood x prior 
evidence 



(10) 



The model's posterior probability function (given the 
data) is explored using MCMC sampling. 

The likelihood of the Schechter para meters is deter- 



mine d using the standard STY method (| Bandage et al 



1979f) . derived from the luminosity function formalism. 
The mass function prior probability is assumed to be a 
top-hat (i.e. weak) prior in the Schechter parameters, with 
boundaries sufficiently wide as to have negligible effects. 
The only adopted exceptions are that M* and $* are as- 
sumed to be non-negative (as we are dealing with physi- 
cal stellar masses), while a < 0. The latter assumption is 
made in order to prevent $(M) deflections in the low mass 
domain, where no constraints are available for SWIRE IR- 
peakers. 

The hyper-parameters are constrained solely by the 
prior knowledge of the stellar mass probability distribu- 
tion function (PDF, provided by the SED fitting); there- 
fore their marginalization automatically accounts for the 
mass uncertainties. The prior probability is thus computed 
as 



Prior ( 



(11) 



We use a Metropolis-Hastings algorithm as our MCMC 
sampler. In order for this sampler to converge efficiently, 
we use a combined strategy for generating proposal steps. 
For the Schechter parameters, we draw new steps from 
univariate Gaussian distributions, whose widths were 
given by dummy MCMC runs. The galaxy mass sampling 



for each galaxy is made by drawing randomly an input- 
sample from the set for each galaxy. We use rejection sam- 
pling, based on the Ay^ of each input sample, relative to 
the best fit ( von Neumann| [l951i) to improve the efficiency 
of this part of the sampling. We note that this second part 
of the sampling is independent of the current position in 
parameter space; this makes the Metropolis-Hastings ac- 
ceptance criterion easy to assess. 

Because each sample is quick to evaluate, we are able 
to calculate a very large number (10^) of steps in one hour 
of 3 GHz CPU time. 



6. 7. Results of parametric analysis 

The results of the Schechter parametric analysis of the 
observed comoving number density of 5.8/Ltm-peakers are 
shown in Fig. [HI based on the Padova-94 library. 

By considering the Schechter sub-space using a 
Bayesian approach, we automatically marginalize over the 
hyper-parameters, producing the posterior probability dis- 
tributions for the mass function parameters. 

The top left panel in Fig. [18] reports the 1, 2, 3 a 
contours in the (M*, a) space, as computed to include 
the 68.3, 95.5 and 99.7 % of the total volume occupied by 
the explored samples. 

The median values obtained for the Schechter charac- 
teristic mass M* are 1.66 x 10" and 1.32 x 10" [ hftf 
Mq], in case the evolutionary correction has or has not 
been considered. The 3a uncertainty range is ~ 1 dex in 
both cases. 

As mentioned above, we are now probing only the very- 
massive tail of the galaxy stellar mass function, hence a 
very loose constraint can be set on the a parameter. We 
find a median value of —0.30, which in fac t is reasonably 
similar to the value a — —0.35 obtained bv lFontana et al ' 



(|2006h on a GOODS-MUSIC sample ranging from - 10^ 
to 4 X 10"'^^ M0. If no evolutionary correction is applied, 
the MCMC space exploration highlights the presence of 
a secondary solution at a = —0.16 (see top left panel in 
Fig. [T8|l . The range of suitable models extends between 
a = -0.12 and -0.40. 

We therefore decided to fix the value of a to —0.35 and 
proceed to explore the (M*, <&*) space (upper right panel 
in Fig. fTS)) . The parameter <&* has a most probable median 
value of 0.00022 [ h^Q Mpc^"^], with a la range smaller 
than 0.1 dex, but doubles if no evolutionary correction is 
considered. Finally, M* is similarly distributed as in the 
former case with a as a free parameter. 

The bottom panels in Fig. [T8| overlay the Schechter re- 
sults on the actual observed 5.8/xm-peaker comoving num- 
ber density, with and without the correction for evolution- 
ary effects in the z = 2 — 3 redshift range. The shaded areas 
belong to the models obtained varying the mass function 
parameters within the la (dark) and 3a (light) ranges, as 
derived with our Bayesian approach. The difference be- 
tween the two cases is very significant, particularly as far 
as M* and $* are concerned. 
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Fig. 18. Results of Schechter parametric analysis of the observed comoving number density of 5.8/im-peakers (Padova- 
94 case). The two top panels report the distribution of models in the M* , a, $* space explored by MCMC sampling, 
in the case of no evolutionary correction. Contours belong to 68.3, 95.5 and 99.7% confidence levels. The bottom plots 
show the possible Schechter mass functions overlaid on the observed comoving number density of 5.8/im-peakers, as 
obtained without (left) and with (right) correction for evolutionary effects. The shaded areas represent the la and 3cr 
ranges obtained in the fit. Error bars include the probability of galaxies to be shifted from one mass bin to another, 
due to uncertainties in the stellar mass estimate. 



The error bars shown in this Figure account also for the 
possible shifts of galaxies from one mass bin to another, 
caused by the uncertainty on stellar mass, as derived from 
SED fitting. The error bars are computed taking into ac- 
count the probability of a galaxy to be shifted, which is 
given by the actual distribution of all ^ 10^ solutions 
in the exploration of the SED parameter space. 

The same MCMC Schechter analysis has been per- 
formed also on the ga l axy m ass function obtained by 
adopting the Maraston ( 2005f) SSP library. The median 



values of the Schechter parameters are very similar to 
those obtained with the Padova-94 library. 

Table [4] summarizes the results, including the median 
and the la, 3cr ranges for t he Schechter par ameters, in 
both the Padova-94 and the Maraston ( 20051 ) cases, and 
with or without evolutionary correction. 



7. Discussion 

We have taken advantage of the wide area offered by 
the SWIRE survey and the extensive multiwavelength 
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Fig. 19. Comparison between the observed z~2 — 3 stellar mass function of 5.8/xm-peakers and literature data, trans- 
formed to a Sal peter IMF (in th e range 0.1 — 10 Mq). Symbol s represent the observed comoving galaxy number den- 
sities from Fontana et al. ( 2006h . Gwvn fc Hartwi ck (200^ and Drorv et al . (2005). The two thick, long dashed curves 
belong to the Schechter STY fit to the 5.8/zm-peakers, as obtained with (lower line) and without (upper line) the evolu - 
tionary correction. The thick, sort-dashed line represents the Sc hechter fit to z = 2 — 3 galaxies in Fontana et al. ( 2006f ). 
The thick dotted line is thc local galaxy stellar mass function bv lCole et al. I (I2OOII V. Thin lin es represents the predictions 
from semi-analytic ([Bower et al.ll2006l iMenci et al.lbool ) and hydro-dynamical ( Nagamine et al. 2005allbh models; the 
solid- shaded area belongs to the Millennium Simulation in the concordance A-CDM cosmogony (jKitzbichler fc White 
200i). 



coverage in the ELAIS-Sl field to identify very massive 
(M > 10^^ M0) galaxies at high redshift. Thus the very 
massive tail of the z = 2 — 3 galaxy stellar mass function 
has been sampled with unprecedented detail. 

7.1. The stellar mass function at z = 2 — 3 

Figure [19] compares the observed comoving number 
of 5.8/im-peake rs to 2: — 2 — 3 data drawn from 
the literature (iFontana et al. 2006t Drorv et al. 



20051 



Gwvn fc Hartwick! 12005^ and transformed to a Salpeter 
IMF (0.1 — 100 Mq). Only Poisson error bars are shown 
here, in order to directly compare our results to literature 
estimates. 

Our data are quite consistent with previous deriva- 
tions of the stellar mass function, at the bright end, but 
deviate in the lower mass bins. This is mainly due to fact 
that no evolutionary correction is usually applied in the 
literature. For comparison, we show the IR-peakers data 



as obtained both with (filled triangles) and without (open 
triangles) this correction. The two thick long-dashed lines 
represent the t wo co rresponding Schechter best fits. The 
Fontana et al. ( 20061 ) estimate (short dashes) is very sim- 



ilar to ours, when no evolution is taken into account. 

SWIRE has the advantage of probing the highest 
masses in an enormous cosmic volume, approaching 10^^ 
Mq, a regime where no previous studies have ever suc- 
ceeded because of the rarity of sources. On the other hand, 
SWIRE is a shallow survey and not enough information is 
available in the low mass regime. As a consequence, below 
Al* the SWIRE IR-peakers Schechter function is poorly 
constrained. 

A significant evolution with respect to the local galaxy 
stellar mass function is confirmed in the highest mass bins 
(M > 5 X 10^^ Mq). This evolution is of the same or- 
der to that derived at lower masses from literature data: 
the number of galaxies at z = 2 — 3 has decreased by a 
factor of ^ 10 at all mass regimes, with respect to the 
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) on the two models and their differences are found in 





median lo- 


3(7 


Menci et all (2006h. 


log(M*) 


ll. 12 11.07 to 11.13 


11.01 to 11.19 


The Nagamine et al. 



$* 0.00048 0.00046 to 0.00052 

a -0.29 -0.36 to -0.17 



0.00032 to 0.00056 
-0.40 to -0.12 



Padova-94 (evo.) 



median 



Icr 



3(7 



log(M*) 11.22 11.18 to 11.24 

$* 0.00022 0.00019 to 0.00025 

a -0.30 -0.34 to -0.18 



11.14 to 11.27 
0.00018 to 0.00031 
-0.37 to -0.12 



Maraston 



median 



1(7 



3(7 



log(M*) 11.10 11.04 to 11.12 

$* 0.00050 0.00046 to 0.00053 

a -0.34 -0.37 to -0.14 



10.97 to 11.15 
0.00035 to 0.00057 
-0.41 to -0.11 



Maraston 



median 



1(7 



3(7 



log(M*) 11.18 11.16 to 11.23 

$* 0.00025 0.00019 to 0.00027 

a -0.34 -0.39 to -0.19 



11.12 to 11.31 
0.00017 to 0.00034 
-0.42 to -0.15 



Table 4. Results of Schechter STY analysis of the comov- 
ing number density of 5.8/.tm-peakers. 



ulations include radiative heating and cooling, super- 
_novae feedback, and standard star formation recipes. The 
two models have been obtained with two different ap- 
=proaches. In Nagamine et al. (2005a) a Smoothed Particle 
.Hydrodynamics (SPH) entropy formulation was adopted, 
-while the other case is based on a Eulerian mesh code with 
a Total Variation Diminishing (TVD) scheme. The SPH 
simulation takes into account also feedback by galactic 
-winds and a multi-component interstellar medium influ- 
encing the star formation process. 

= The agreement of these four models with the data is — 
-unfortunately — rather poor. The two hydro-dynamical 
-simulations systematically overpredict the number density 
of galaxies at z = 2-3, at 10^° < M < 10" M©, as shown 
by comparison to literature data. Moreover the SPH run 
"applies a cutoff at M = 5 x 10"'^"'^ M©, while our data also 
show a non-null number density at these very high masses. 

The SAM predictions are closer to the actual data, 
but merely intersect the observed $(M), under- predicting 
the comoving number density of very massive objects and 
overpredicting that of lower- mass sources. This happens 
_at M ~ 10" M0 for the Menci et al. (2006) model and at 
M ~ 3 X 10" Mq for lBower et all (|2006t ). The latter shows 
a very steep slope above M* , too steep to be consistent 
with the observational evidence. 



current epoch. On the other hand, at lower redshift only a 
weak evolution is detected for m assive galaxies {M > 10^ ^ 



Mr^), at 



Bundv et al 



leas t to z 
20051: 



1-5 (e.g. 



Fontana et al 



Franceschini et al.l l2006t 
20061 among others) , 



while the number of lower mass objects is significantly 
lower at z ~ 1 than the local value. 

The solid-shaded area in Fig. [12] represents the predic- 
tion by the Millennium A-CDM Simulation, and is taken 



from Kitzbichler fc White llooi). The wide area covered 



by the model is meant to represent measurement errors 
in log {Mi,). The results of this simulation are fairly con- 
sistent with our 5.8^m-peaker data, above 2 x 10^^ M©. 
Nevertheless, the slope of their prediction appears to be 
too steep, because the number of 2; = 2 — 3 galaxies be- 
tween M = 1 — 2 X 10^^ M0 is overpredicted. Note also 
that below 10^^ Mq, the galaxy stellar mass function is 
significantly overestimated by the Millennium model, with 
respect to the observed literature data. 

Thin lines and other shaded areas belong to s emi- 
analytic (SAM, iBower et a l."2006': 'Menci et al.ll2006h and 
hydro-dynamical In^ gamine et al.i,2005aK ,b) models. Both 
the two SAMs include AGN feedback on star formation, 
although in two di fferent ways: in the Durham simulation 
( Bower et al ][200i) the AGN feedback is related to the gas 
smooth accretion onto the central black hol e, continuing 
into t he AGN quiescent phase; while in the iMenci et al.l 
(200i) case it is produced by shock waves originated dur- 
ing the short active AGN phase only. Further details 



7.2. The integrated stellar mass density 

A complete view on the contribution of high-redshift mas- 
sive galaxies to the mass budget in the Universe is given 
by their global stellar mass density, as obtained by inte- 
grating the stellar mass function. 

As far as 5.8/im-peakers are concerned, we have de- 
rived this quantity by exploiting the properties of the 
Schechter function. Its integral can be expressed as: 



<^{M)dM 
= M*$* X F ( a + 1 



M, 



inf 



M* 



(12) 



where T{a,x) is the incomplete gamma function, esti- 
mated from a to infinity, and Minf is the lower-mass inte- 
gration cutoff (set by the completeness limit, in our case). 
We integrate the stellar mass function over the mass range 
for which we were able to apply a reliable completeness 
correction, i.e. M > 1.6 x 10^^ Mq, taking into account 
the large uncertainties in the Schechter parameters, due 
to poor sampling in the low mass regime. The resulting 
stellar mass density (accounting for evolutionary correc- 
tion) is yO* = 1.18 X 10^ [hyo M0 Mpc~^], with a 3cr range 
between 6.0 x 10^ and 2.3 x 10^ [hro Mq Mpc'^]. 

In the case of 4.5/im-peakers, instead, we can set only 
a lower limit to the actual stellar mass density built in 
M > 10^^ Mq. By simply integrating the observed data, 
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Fig. 20. Integrated stellar mass density as a function of redshift. The two large symbols represent the values of Pi, 
derived from our data. The filled triangle represents the contribution to the stellar mass density by 5.8/im-peakers 
above 1.6 x 10^^ M0, as obtained by integrating the Schechter solutions of our MCMC fit. The filled circle sets a 
lower limit to the stellar mass density of 4.5/im-peakers above 1.25 x 10^^ M©, obtained by integrating their observed 
number density. Literature data from numerous works are shown (see right panel for references) , integrated down to 
10^ M0. Thin lines belo ng to different models , whi le the thick solid l ine simply connects the obtained for galaxies 
heavier than 10" Mq bv lFontana et all (|2006f) and lCole et"all (|200l[ ). All data and models have been transformed to 
a Salpeter IMF, extended to the mass range 0.1 — 100 Mq. 



we obtain > 6.55x10^ [hyo Mq Mpc'^]. Table[5]reports 
the stellar mass densities thus obtained. 





4.5/im-peak 


5.8^m-peak 




P* [h-70 


Mq Mpc-3] 


Observed {Va, no evo) 


> 6.55 X 10" 


1.45 ±0.3 X 10' 


Schechter STY fit (no evo) 




1.56 X 10^ 


Schechter 3cr (no evo) 




0.5 - 2.9 X lO'^ 


Observed {Va, evo) 




1.23 ±0.3 X 10'' 


Schechter STY fit (evo) 




1.18 X 10^ 


Schechter 3(7 (evo) 




0.6 - 2.3 X lO'^ 



Table 5. Stellar mass density for 4.5/im-peakers {z = 1 — 

2, M > 1.25 X 10" Mq) and 5.8/im-peakers (z = 2 - 

3, M > 1.6 X 10" Mq), as obtained with the 14 and 
parametric analyses (Padova-94 library). 



Figure [20] compares the based on IR-peakers to 



2005 



2005 



2no4r I Franceschini et al 
Gwvn fc Hartwicd l2o'o5 



a col l ection of data from the literature (IRudnick et al 
2006L l2003l iFontana et allliooel l2004l. l2003HDrorv eral 



I2OO6I: iBundv et al 



Caputi et al.l 



2005[ iDickinson et all l2003bl: iBrinchmann fc Ellid boOC : 



2006 . 



Cole et al.l 120011 ). all transformed to a Salpeter IMF 
extended between 0.1 and 100 Mq. This graph highlights 
a dramatic scatter in the current estimate of the stellar 
mass density, with significant discrepancies between the 
various authors. 

While the literature data were obtained by integrat- 
ing the mass function down to 10* Mq, the IR-peaker 
data points belong to sources above ~ 10^"'^ Mq only. The 
solid thick line simply connects the p^, estimate fo r galax - 
les more massive than 10" Mq bv lFontana et all (|2006l ). 



Despite the large uncertainty on the actual value, due 
to degeneracies in the Schechter fit, our measured_£4_Jor 



the 5. 8/jm-peakers is fully consistent with iFontana et al 
( 20061 ) data. On average between 30% and 50% of the 
total stellar mass in galaxies at z = 2 — 3 is stored in 
our population of massive (M > 1.6 x 10^^ Mq) 5.8pm- 
peakers. Although v ery unlikely (on the basis of data by 
IFontana et al. 20061 ). this value could in fact grow to a 
higher fraction if we account for the uncertainties and/or 
compare to the lowest estimate of the overall stellar mass 
density from the literature. 

The thin lines in Fig. [201 represent the evolution of the 
global stellar mass density, as predicted by the same semi- 
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ana lytical and hydro-dynamical model s shown in Fig 
T^l (iBower et al.ll2006l : iMenci et all 120061 : iNagamine et al, 
2005allhr i. For all models the evolution of the stellar mass 



density becomes steeper at z > 1.5 — 2.0. The values of 
predicted by different models gradually diverge at redshift 
z > 1, therefore the observational data could in principle 
constrain the actual scenario, but the current wide scatter 
and large error bars in the stellar mass density makes it 
very hard to disentangle the various models. However, we 
provide the first tight constraint at the largest masses, for 
z = 2-3. 

The top panel of Fig. [21] shows the ratio of the stellar 
mass density at a given redshift and in the local Universe, 
for different mass bins. Crosses represent the data ob- 
tained by integrating the Schechter fit to the mass function 
described by Fontana et al. (2006,, see also their Fig. 9). 
We have chosen these authors, because their analysis ex- 
tends to the redshift bin covered by 5.8/im-peakers and 
because they provide Schechter parameters in all redshift 
bins. The dotted, dashed and solid lines represent three 
different mass bins. Triangles represent 5.8/j,m peakers, 
and are obtained by integrating our Schechter fit in the 
10" < M < 5 X 10" and 5 x 10" <M< lO^^ M© mass 
bins. 

Despite the poor statistics and the high masses and 
the large "noise" p reviously pointed out in Fig. [20l 
iFontana et al. ( 2006f ) data clearly show that very massive 
galaxies evolve more rapidly than lower mass objects and 
reach the local density at earlier epochs. It is also inter- 
esting to note that 5.8/im-peakers sample a redshift range 
where /5*//5* (0) assumes similar values in all the three 
mass bins considered. This effect explains why the ratio 
of the stellar mass function at z ~ 2.5 to the local one is 
relatively flat (bottom panel of Fig. [2T|) . The latter plot 
shows that spans values within a factor of 3 — 4 in 

the whole mass range between 5 x 10^ and 10^^ M©. 

At least 50-60% of the present-day stellar mass in very 
massive systems seems to have assembled between z ~ 2.5 
and z = 1.5, i.e. when the Universe was between ^2.5 and 
~ 4.5 Gyr old. This process was roughly complete at z ~ 
1, when more than 80% of massive galaxies were already 
in place. SWIRE 5.8/im-peakers (z = 2 — 3) represent less 
than ~10% of the stellar mass in the local population of 
massive galaxies (M = 10" - 10^^ M©). 

According to theoretical models, the formation of 
galaxies and large scale structure occurs in the frame of 
some variant of 'bias ed' hierarchical buildup within a A- 
CDM cosmology (e.g . ICole et al.llioool: iHatton et aLll2003t 



Granato et al.l |200j) . In these scenarios the most mas- 



sive objects (e.g. > several 10^^ M©) are predicted 
to assemble earlier, more quickly and in richer environ- 
ments than less massiv e ones (e.g. ISomerville et al. 2001 
Nagamine et al.ll2005bl) . 



It is worth noting that no information is currently 
available on the environment hosting these very massive 
galaxies at high-redshift. The study of the clustering prop- 
erties of mid-IR powerful-emitting IR-peakers has shown 
that their spatial correlation length is rg = 14.40 ± 1.99 
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Fig. 21. Top panel: ratio of the stellar mass density at 
a given redshift and in the local Universe, as a function 
of redshift and for different mass bins. D atapoints with 



crosse s were computed by integrating the iFontana et al 



( 20061 ) Schechter fit to the stellar mass function. Bottom 
panel: ratio of the mass function at z = 2 — 3 to the lo- 
cal one, as derived for 5.8/im-peak ers (triangles) and from 
the literature (other symbols). The I Cole et al.l (|200l[ ) local 
mass function was adopted. 



[h^i Mpc] at z = 2 - 3 (jFarrah et al.ll2006l ). This value is 
consistent with a populations residing in M ~ 10"'^^ M© 
dark matter haloes, thus they are excellent candidates for 
being part of protoclusters at high redshifts. If this is the 
case, and if less active IR-peakers (i.e. with lower mid-IR 
luminosities or not detected by MIPS at all) showed anal- 
ogous properties, they could trace quiescent stages of simi- 
lar systems. Broadly speaking, the studied 5.8/im-peakers, 
sampling the very massive tail of the stellar mass function 
at z = 2 — 3, might represent a different class with respect 
to the field general population. They could, in fact, be 
at the center of large dark matter haloes, their compan- 
ions being too faint to be detected by SWIRE. Different 
environmental conditions would rule their evolution and 
partially explain the mass "downsizing" effect exemplified 
by a faster high-z evolution of massive systems. 



28 



Berta S., et al.: IR-peakers mass function 



8. Summary and conclusions 

Sampling the very massive tail of the stellar mass func- 
tion at high redshift and estimating its contribution to 
the global stellar mass density is a critical task in mod- 
ern cosmology, motivated by the recent evidence that a 
"downsizing" effect exists in the evolution of stellar mass 
across cosmic time. 

We have exploited SWIRE/Spitzer and ancillary data 
in the ELAIS-Sl area to perform a systematic search for 
high redshift (z > 1.0) massive {M > 10^^ Mq) galaxies. 
High redshift systems have been isolated by identifying 
the 1.6/im restframe stellar peak shifted to IRAC wave- 
lengths (3.6—8.0 /im). The availability of near-IR (J and 
Ks band) data allowed us to avoid low-redshift interlopers 
and selection aliasing due to bright 3.3 /im PAH features 
at z ~ 0.4. A total of 203 5.8/im-peakers and 123 4.5/j,m- 
peakers have been identified over one square degree in the 
ELAIS-Sl field. 

We have performed an extensive SED analysis, based 
on mixed stellar population synthesis, focused on deriving 
the stellar masses of the selected sample. The advantage of 
near-IR and mid-IR constraints, as well as the dependence 
of results on the choice of the IMF and SSP library have 
been explored in detail. The main results are: 



— because of the shallow 5.8/im flux cut adopted, the 
SWIRE IR-peaker sample consists of very massive 
galaxies, the majority of sources having > 10^^ 
Mq . Objects in the range z = 1 — 2 peak at Af^ ~ 10^^ 
Mq, while the distribution of 5.8/j,m-peakers (z = 2—3) 
is centered at M^, ~ 2 x 10^"'^ Mq. Typical uncertain- 
ties in stellar mass estimate (due to degeneracies in the 
SFH space) range between 0.1 and 0.3 dex, depending 
on multiwavelength coverage. The emission of ~30% of 
the sources turns out to be dominated by stars older 
than 1 Gyr, and also in the majority of the remain- 
ing cases old stellar populations do contribute to the 
observed SEDs. 

— the availability of mid-IR data provides a valuable con- 
straint on the recent star formation history of individ- 
ual galaxies. If no 24/xm flux (nor upper limit) were 
available, the resulting stellar masses could be under- 
estimated and the spread in mass would be wider. 

— despite being very useful in the selection process, near- 
IR {J and Kg) data turned out to be effective in con- 
straining the D4000 break only in ~15% of cases, and 
preferentially fo r 4.5/xm-p e akers. 

— the choi c e of a ChabrieJ ( 2003al ) IMF, instead of a 
ISalpete one, leads to systematically lower stel- 
lar masses, resulting in a rigid shift of the stellar mass 
distribution by ^ 0.3 dex. 

— when includi ng thermally-pu lsing AGB stars in the 
SSP library ( Marastonl 200^, the fraction of objects 
dominated by old stellar populations increases, but the 
overall stellar mass distribution does not change signif- 
icantly, because of the different M^/L ratios of SSPs. 



The stellar mass estimates have been used to compute the 
comoving number density of galaxies as a function of stel- 
lar mass (i.e. the observed stellar mass func tion), adopting, 
the ac cessible volume formalism. Following iFontana et al.l 
(j2004l ). we have corrected the samples for mass incom- 



pleteness and recovered the mass function down to 1.25 x 
10^^ and 1.6 x 10^^ Mq for 4.5/Ltm- and 5.8/.tm-peakers 
respectively. 

Unfortunately, the selection, based on a 5.8fim flux cut 
turned out to be only partially sensitive to 4.5/im peakers, 
therefore only a lower limit on the mass function of z = 
1 — 2 sub-sample could be set. 

The observed stellar mass function of 5.8/im-peakers 
was reproduced with a parametric function, using the STY 
ISandage et al. 19791 ) approach. The uncertainties in the 
stellar masses of individual sources were automatically in- 
cluded in the analysis by using a Bayesian formalism, and 
the best fit was obtained with a MCMC sampling of the 
parameter space. The stellar mass function was finally in- 
tegrated to derive the stellar mass density locked in 5.8/im- 
peakers with M > 1.6 X 10" Mq, at z = 2-3. The resuhs 
of the mass function analysis are: 

— the wide area surveyed by SWIRE allows the very mas- 
sive tail of the stellar mass function to be probed with 
unprecedented detail at z = 2 — 3, extending previous 
analyses up to = 7 x 10" [hf^f Mq]. 

— at M < 5 X 10^-'^ M0, a significant intrinsic evolution 
has been detected across the redshift range z = 2 — 3, 
strongly dependent on mass. The dependence or the 
5.8//m-peakers number density on (1 + z) has powers 

of 0.4 and 0.65 for M - 2 x 10" and 4 x 

10^^ M0 respectively. In the highest mass bins (M > 
5 x 10^^ Mq) the number density of IR-peakers keeps 
nearly constant. 

— comparison to literature data for the z = 2 — 3 mass 
function shows an overall agreement of the 5.8/im- 
peaker comoving number density to that of the K20, 
MUNICS, and GOODS-MUSIC surveys, in the higher 
mass regime. At lower masses a significant evolution- 
ary correction should be applied, when the mass func- 
tion is averaged over a wide redshift range. 

— a significant evolution of the stellar mass function of 
M ^ 10^^ Mq galaxies with respect to the local esti- 
mate was detected: $(z = 2 - 3) < 0.1 x <i>(z = 0). 
Combining 5.8/xm-p eakers and literature data (e.g. 
Fontana et al.|[2006f) . this implies that the bulk of mas- 



sive galaxies was not yet in place by the time the 
Universe was ^ 3 Gyr old, but must have been as- 
sembled in the following ~ 1.5 Gyr of evolution. 

— current hydro-dynamical models significantly overesti- 
mate the number density of massive galaxies, while the 
semi-analytic approach underestimates it. 

— since SWIRE 5.8/im-peakers sample only the very mas- 
sive tail of the mass function, the Schechter slope a 
cannot be constrained. Despite its very low signifi- 
cance, the best fit value a — —0.30 is similar to 
that found in the literature. The other best fit pa- 
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rameters are: M* = 1.66 x 10" [hf^^ Mq] ±1 dex; 
$*= 0.00022t^;S[h?oMpc-3]. 

— the integrated stellar mass density of 5.8/xm peakers is 

= 1.18 X 10'^ [li7o Mq Mpc~^], with a 3a range of 
±0.3 dex. Only a lower limit could be set for 4.5/um- 
peak galaxies: > 6.55 x 10*^ [hro Mq Mpc"^]. 

— on average SWIRE massive 5.8/xm-peakers provide 
30—50% of the total stellar mass density in galaxies 
at z = 2 — 3. 

— 5.8/im-peakers provide less than 10% of the stellar 

mass locked in massive galaxies (M = lO"^-"^ — 10^^ Mq) 
in the local Universe. 

The analysis carried out on SWIRE massive galaxies at 
z > 1.5 over one square degree, highlighted the comple- 
mentarity of wide-shallow and deep pencil-beam surveys. 

On one hand, sampling the faint end of the luminos- 
ity function, i.e. the low-mass end of the mass function, is 
needed in order to constrain the shape of the mass func- 
tion and the total stellar mass density in galaxies at high 
redshift. 

On the other hand, very massive galaxies are rare ob- 
jects on the sky, with a number density 1.2 x lO"** [h^Q 
Mpc"'^] and it is necessary to explore large volumes of 
Universe in order to fully characterize them. The natu- 
ral extension of this analysis is to build the stellar mass 
function of high- 2; galaxies over the whole SWIRE 49 deg^ 
area, taking advantage of what we have learned thanks to 
the full multi-wavelength coverage in ELAIS-Sl. 

At the time being, not much information is known 
about the environment hosting these massive high- 2: ob- 
jects. Further analyses of this popiilation should be carried 
out probing also the environmental frame they belong to, 
in order to correctly interpret their role in the "downsiz- 
ing" scenario. 
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